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THE LARGE-SCALE CRUMPLING OF THIN 
CYLINDRICAL COLUMNS 
By SIR ALFRED PUGSLEY and M. MACAULAY 
(The University, Bristol) 
{Received 11 December 1958] 


SUMMARY 


When a very thin metal tube of cylindrical section is compressed between parallel 
platens, its walls tend to buckle in and out to form a diamond pattern of deforma- 
tion around the tube. This paper considers the subsequent behaviour of such tubes 
when the compression is continued to cause large-scale crumpling of the tube walls. 
The nature and mode of this crumpling is examined in the light of experimental 
results and, guided by an approximate theory for an idealized case, an empirical 
expression for the load required to effect the crumpling is obtained, 


1. Tue elastic buckling of very thin tubes subject to axial compression 
has been the subject of extensive study (1, 2), both theoretical and 
experimental, during the first half of this century. As a result, it has 
become well known that for such cylindrical columns: 

(a) The critical buckling load, as determined by experiment, is much 
below (roughly one-third of) that given by the classical theory based on 
small deflexions. 


(6) The local buckling of the tube walls, as observed by experiment, 
does not conform in size or pattern with the results of classical theory. 
In general, the walls buckle into a diamond pattern,+ each diamond being 
approximately square with its diagonals in the longitudinal and circum- 
ferential directions. The diamonds so formed tend to be appreciably 
larger than suggested by classical theory and to favour inward rather 
than outward deformation. 


The departures from classical theory represented by (a) and (b) are 
now thought to be due partly to initial geometrical irregularities in the 
curvature of the tube walls, and partly—perhaps largely—to the non- 
linear support that longitudinal elements or strips of the tube walls receive 
from the circular circumferential strips. A circular ring provides decreasing 
resistance to radial compression as tho inward deflexion proceeds. On 
this basis, a physical explanation of the preferential inward deformation 
of a thin cylindrical column is forthcoming, together with a qualitative 
explanation of the lowering of the critical buckling load and the increase 
in the size of the individual ‘diamonds’. 

2. Most of the experimental work done to date has concentrated on 


+ The symmetrical, concertina-like buckling characteristic of rather thicker tubes is not 
dealt with here. 

(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960} 

5092.49 B 
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the features (a) and (b) above, and as a result few compression tests 
on thin tubes have been carried to large axial strains. In recent years, 
however, the authors have been concerned} with studying the energy 
absorbed by thin tubes when they are compressed by amounts sensible 
in relation to their radii and/or lengths. When this is done, following upon 
initial buckling into the well-known diamond pattern, these buckles first 
become accentuated and then tend (at least over a length of the tube) to 


Fic. 1 (b) 


increase in size by ‘swallowing’ adjoining buckles, usually until a final 
stage is reached when there are only three or four large diamonds around 
any circumference of the tube. Further compression then ‘folds’ these 
large diamonds to a closed pattern, so shortening the column in the region 
of the folds very considerably. This large-scale crumpling process is 
illustrated by the pictures given in Fig. 1. 


The purpose of the present paper is to obtain a physical understanding 
of the crumpling process and thence, guided by a theoretical analysis of 
a simplified version of this process, an empirical expression for the work 
done in crumpling a thin tube under axial compression; and, incidentally, 
also an estimate of the average resistance to crumpling. It is thought that 
the bulk of the work done must arise from plastic deformation of the metal 


+ With the support of the Aluminium Development Association and the British Transport 
Commission. 
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of the tube walls; no account, therefore, will be taken of any residual elastic 
strain energy stored in the tube. 


3. We have first to consider the geometry of the crumpling process,t and 
so if possible establish an idealized process embodying the main features 
of the practical deformations illustrated by Fig. t. When the large-scale 
diamonds are established, and before any further folding or crumpling 
has occurred, the tube, originally circular in section, will appear essentially 


Fic. 2 (a) Fia. 2 (b) 


in one or other of the forms shown in Fig. 2, accordingly as there are three 
or four diamonds (n = 3 or 4) around its circumference. At this stage, 
a little plastic work will have been stored along the rather sharp edges 
and diagonals of the diamonds and the tube will have shortened some- 
what (19 per cent for n = 3; 8 per cent for n = 4). At this stage too, 
the plan form of the tube will appear as in Fig. 3, and if no appreciable 
stretching of the tube material has occurred it is evident that if d is the 
length of a circumferential diagonal, then 

nd = 2nR, (1) 
where # is the radius of the undeformed tube. 

With further compression the large ‘square’ diamonds collapse on them- 
selves by folding across their circumferential diagonals and along their 
sides. In doing so they tend to maintain their initial plan form, so that, 
for example, the tube of Fig. 2 (a) and Fig. 3 (a) appears, when its diamonds 


+ Since this paper went to print, the authors’ attention has been drawn to an interesting 
early contribution to this subject by A Mallock (3). 
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are folded flat, with a plan geometry still approximating to that of Fig. 
3(a). The great resistance of the metal to circumferential extension has, 
in effect, caused the crumpling process to be executed as though the metal 
walls were guided axially down rods at each of the corners of the hexagonal 
pattern of Fig. 3(a). To a first order, therefore, the plan form of the tube 
when the large diamonds are first formed is also that of the tube when 
fully crumpled. 

In this process, if the folding were per- 
fect, the actual area of the metal walls 
would necessarily have been greatly re- 
duced. Such an impossibly large reduc- 
tion of area (it would be down to 59 per 
cent of the original area for n = 3 and 
41 per cent for n = 4) is in practice 
avoided by the fact that the diamonds do 
not fold really flat along their edges and 
circumferential diagonals. A section at 
right angles to an edge or diagonal shows 
a finite region of circular bending and not 
a sudden localized fold; this is, of course, 
consistent with experience when one tries 
to bend any strip of metal locally—it 
develops a finite curvature in preference 
to a sudden kink. As a result, an appre- 
ciable part of the original wall area has, 
when a tube is crumpled nearly flat, gone 
into the roll-like edges formed, as well as 
into the triangular shaped flat areas. 

It is necessary to notice also that this restricted folding process, coupled 
with curvature at the edges, cannot occur if the metal remains continuous 
without considerable shear deformation of the metal in its plane, parti- 
cularly in the region of the corners of each of the flat triangles formed in 
plan. On the simplest view, each large diamond, originally of square shape, 
has, by the crumpling process, become a folded rhomboid as though 
sheared in the manner indicated in Fig. 4. However, on account of the 
way this process is restricted circumferentially, as already explained, and 
because it is complicated by the resulting formation of ‘rolled’ edges to 
the triangles formed, this simple shear view of the deformation is very 
rough indeed and the large-scale plastic deformation in shear, as indicated 
in Fig. 5, tends to be concentrated in the corner regions. In these regions 
also it is clear that some plastic extension, as well as shear, must occur. 
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4. Whatever the number of folded diamonds that are finally produced 
around the tube circumference, each diamond will have started as a square 
one and ended as a rhomboid folded about its circumferential diagonal. 
The effective length of this diagonal, initially given by equation (1) above, 
will remain approximately the same throughout the crumpling process; 
but the angle between it and any side, initially about 45°, will fall to 30° 
when x 3 and 22)° when n 4, as illustrated in Fig. 3. 


On this simple geometrical basis we can now assess the energy absorbed by 
the deformation of the diamond pattern to the folded rhomboid pattern. 

First as to the work done in plastic bending along the lines of the final 
folds. Coupled with the folding along a diagonal of length d there will 
also be folding along each side s of the diamond. If we consider a length 
of the undeformed tube equal to 4d, corresponding to half the height of 
the original square diamonds, there will be foldings along two sides to 
be accounted for per diamond in addition to one diagonal fold, making 
a fold of length (d+-2s) in all. Now the work done per unit length of fold 
will be given by the product of the full plastic moment required to bend 
the wall plating (per unit length) freely and of the angle (7 radians) through 
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which it is turned. On this basis, if p is the yield stress for the wall 
material the above plastic moment is }pt? and the total plastic work done 
in the folding process around the tube is given by 

= 4(mn)pt?(d+-2s) = kh, npt*d, (2) 
where /, is a coefticient depending on the relation between s and d. For 
n = 3, 8 = 0-58d, and k, = 1-70; for n = 4, s = 0-54d, and k, = 1-64. 

If we now consider the deformation of the shape of the diamonds from 
square to rhomboid this, as indicated in Fig. 4, can be regarded as a plastic 
shear distortion under forces along the sides of each diamond. If q is the 
average shear stress concerned, the force along any side will be gts and 
the work done by this action, for an angle change of @ as between square 
and rhomboid, will be, for the diamonds in a length of tube equal to 4d, 

U, nOqts* k, ngtd®, (3) 
where /, is a coefficient depending on the relation between s and d. For 
n 3, k, = 0-18; for n 4, k, = 0°23. 

If we now write P as the average compressive force acting on this tube 
during the crumpling process, we can equate the internal and external 
work thus: }Pd = U,+U, = ky nptd+k, ngtd?, (4) 
and so derive P = 2k, npt?+-2k, nqtd. (5) 
In equation (5), since k, and k, vary only slowly with n, the first term will 
be least when n is least; and since nd is (from (1)) constant, the second 
term hardly varies with n, though it is least when k, is least, corresponding 
to n being small. But m must be an integer not less than 3; hence P will 
be least when x = 3. With the appropriate values for k, and k, in (5), 
this becomes P = 10-2pt?-+ 1-08qtd. (6) 
It will be convenient to compare this load P with the end load P, that 


would cause simple yielding of the tube in direct compression. With the 


same notation P, ia 2x Rtp, (7) 


and the ratio of these two loads is thus given by 
P 
—=(C, ‘+6, 8 
where C, = 1-62, C= 0-36 2. 
P 
It remains to consider the value of g/p. This ratio of yield stresses, for 
most metals, is of the order of 0-5; but it may well be, because of the 
distribution of shear noted in section 3, that q/p is effectively less than }. 
Taking g/p = } as a tentative estimate, 
P 


t 
-= 16—+40°12. 9 
P. pte! (9) 
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Before proceeding to consider (8) and (9) in the light of experimental 
results, we have first to recall that they are based on n = 3 and to con- 
sider how they would be affected were n larger than 3. For n = 4, for 
example, k, and k, have already been evaluated, and (5) and (7) give at 


once 
P t 
— as 9) —-1- 


-- 0-15 10 
R a" C (70) 
which is substantially greater than (9) for any value of t/R. 


34 


Load (tens). 











+ 


2 


Shortening Cin) 


Fic. 6 


5. When a thin-walled tube is slowly crushed between parallel platens, 
it crumples progressively in the manner illustrated by the photographs 
of Fig. 1. At the same time its resistance varies with the amount of 
shortening of the tube in the manner illustrated by the curve of Fig. 6. 
In order to compare experimental results with expressions such as (8) or 
(9), we have first to integrate the area under the curve of Fig. 6 and thence, 
by dividing by the total shortening of the tube, derive an average value 
for the resistance of the tube to crumpling. This average resistance, based 
on the work done in shortening the tube, will be a direct experimental 
measure of P in our equations. 

This has now been done for a number of tubes with various values of 
t/R and with materials (stainless steel and soft aluminium) for which the 
yield stress p (taken as the 0-2 per cent proof stress) has been measured. 
As a result the ratio P/P, so determined experimentally, has been plotted 
against ¢/R in Fig. 7. Two curves are shown on the diagram against the 
background of the experimental points thus obtained. One corresponds 
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o (9) directly and the other corresponds to (8) with C, and C, chosen to 
give a close fit to the experimental points. It will be seen that (9) itself 
gives a rough approximation to the experimental results and that (8), 
rewritten as P 


P, 


t 
= §— 0-1 
$3418, (11) 


fits very well. 


1-O-4 


© = experimental results 


P/Py = 5 t/R+0.13 


PLP = 16 t/R +0-l2 











6. The foregoing simple theory is given as a first approach to the 
problem of determining the energy absorbed in the crumpling of a very 
thin circular tube in compression. It has led to a form (equation (8)) 
capable of fitting our experimental results, but the theoretical coefficients 
(in equations (9) and (10)), though of the right order, clearly require to 
be replaced by others established experimentally. It is of interest to note 
that in (9) the coefficient of t/R appears more inaccurate than the constant 
term. This may be due partly to theoretical error, arising, for example, 
from the extra work done in folding at the diamond corners, but partly 
also to the fact that the en:pirical value in (11) depends largely upon 
experimental results for high t/R values. Here the mode of buckling is 
approaching a transition to the symmetrical concertina-like form typical 
of rather thicker tubes, for which it is known (from work by the second 
author, for publication elsewhere) that P/P, tends to vary with ,/(t/R) rather 
than ¢t/R. Further work bearing on these points is in hand. 
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AN APPROXIMATE ANALYSIS OF THE COLLAPSE OF 
THIN CYLINDRICAL SHELLS UNDER AXIAL LOADING 


By J. M. ALEXANDER 
(Imperial College of Science and Technology) 


[Received 29 January 1959] 


SUMMARY 
An approximate theory for the process is derived, leading to a solution of the type 
P = Ct'*vD, where P is the collapse load, ¢ the shell thickness, D the shell diameter, 
and C a constant for any given material. Good agreement is exhibited between this 
relationship and experimental results. 


1. Introduction 


Ix the design of nuclear reactors having vertical fuel channels, it is 
necessary to guard against the accidental dropping of components, either 
in these fuel channels, or in other vertical channels containing, for example, 
control rods. To absorb the energy consequent on dropping some com- 
ponent, it is usual to provide energy-absorbing devices. One such device 
takes the form of a thin cylindrical shell arranged to buckle lengthwise 
when hit by the dropped component. 

Considerations of damage to the components set limits on the allowable 
loads which can be sustained during this absorption of energy. It is there- 
fore necessary to determine a relationship which will allow prediction of 
the dimensions of a shell to give the required resisting load. An empirical 
relationship can be determined by experiment—the following note gives 


an approximate theory for the process, to provide a guide for the choice 
of the empirical formula. 


2. Theory 


A simple mode of collapse of the tube is assumed, and the work required 
to achieve this mode determined. The mean collapse load is determined 
by equating this work to the mean collapse load multiplied by the distance 
through which the load operates, 

There are two actual modes of collapse, one in which the tube forms 
symmetrical convolutions so that it takes up the appearance of a bellows, 
the other in which both transverse and longitudinal waves are formed. 
This latter mode is the general mode of collapse for thin tubes, but is too 
complicated to use for the present purpose. A simplified version of the 
‘bellows’ type of collapse mode is adopted, it being assumed that the shell 

[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960] 
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collapses in the form of a ‘concertina’ with straight-sided convolutions. 
One of these convolutions is illustrated diagrammatically in Fig. 1, the 
symbols used for the various dimensions being also shown on that figure. 


Pp 


ll et te 





























a ai al a ae ae ke ea a ar ae a a 


p 
Fic. 1. Assumed collapse mode. 


The work done in deforming the metal into one such convolution can 
be split into two parts, namely that required for bending at the circular 
‘joints’, and that required for stretching the metal between the joints. 

To simplify the analysis, elastic strains and work hardening of the 
material are neglected, i.e. a ‘plastic-rigid’ material is assumed. 

During an increment dé of the angle 6, the increment of work done at the 
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three joints shown is 

dW, = 4M d@n(D-+-hsin 8), (1) 
where M is the collapse moment per unit circumferential length of the 
joint. 

For a narrow beam M would be equal to Yt?/4, Y being the yield stress 
in uniaxial tension or compression (assumed equal). In this case the 
‘beam’ is essentially very wide and may be regarded as being deformed 
under substantially plane strain conditions during the incremental change 
dé. Under these circumstances, if the material obeys the von Mises 
criterion of yielding the direct stress will be raised to (2/¥3)Y, so that the 
collapse moment M will become (2/¥3)(Y?t?/4). Thus 

2a P 
dW, = 3 Yt? d0(D-+-hsin 6). (2) 

The mean strain in extending the metal between the joints, during the 

incremental change d@, is 
7 D+hsin(@+d6)|—7[D+-hsin6] _hd@cos@ 
a(D+hAsin 8) ~ D+hsind 

The stress in these fibres will be equal to Y, the yield stress of the metal, 
so that the increment of work done in stretching them will be 

dW, = 14400089 |. hsin6). 2ht = 2nYh%d0 0080. (4) 
“ D+hsin@ 
The total work done in collapsing one convolution, i.e. for @ increasing 
from 0 to 90°, is therefore 


(3) 





us 


W= | (dW,+am,) | E YD +hsin8)+2aYh* e036] dé. (5) 
a . Ne 
0 


This must be equal to the mean axial load P multiplied by its total 
displacement 2h (neglecting the thickness of the metal). 
Hence, integrating equation (5), 
P : mt? (7D 
rs ‘| 2h” 


h is now determined, by minimizing this expression, to give 


y |aht; 


/ TT , 
h = Jes) = 0-953 ./ (Dt). 


It is interesting to compare this value of h with that derived from the 
elastic analysis for the buckling of thin cylindrical shells,+ namely 


7 (Dt) 


h= sm nT 
2{3(1—v*)} 


(8) 


t+ Timoshenko, 8., Theory of Elastic Stability, p. 441. 
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where v = Poisson’s ratio. Substituting v = 0-25 gives 
h = 1-213,/(Dt). (9) 
Substituting h = k,/(Dt) in equation (6), 


Y= \2v3k 

An alternative mode of collapse is to assume that the convolutions are 
formed internaily instead of externally as shown in Fig. 1. If the analysis 
is repeated for this case there results the following variation of equation 
<= (—_ + wk |e 5vD —= 

Y (; vV3k - ‘ 


Bearing in mind the approximate nature of this analysis, and also that 
the true deformation mode lies somewhere between these two cases, it 
seems reasonable to adopt the mean value between equations (10) and 


(11), i.e. Pp ( 


Dp 2 2 
r= ( 7 +rk)eD 47. (10) 


{2 


5 (11) 


y= 


a +ak|tsvD, (12) 
2V3k ° 


Substituting k = 0-953 from equation (7) gives 
P 
- — 5-991 VD, 


and substituting k = 1-213 from equation (9) gives 


2 
—_= 6-16th)s dD. 


Thus, the final solution becomes of the form 


P= KY#4yD, (15) 
where A = 6-08, Y = yield strength of the material (assumed perfectly 
plastic), ¢ = thickness of the cylindrical shell, D = mean diameter. 


3. Comparison with experiment 

Table 1 gives data relating to experiments carried out to determine the 
mean collapse loads for mild steel tubes. 

The data are shown plotted in Fig. 2, on which the ordinates are P and 
the abscissae t’°.D. The mean straight line drawn through these points 
and the origin has a slope AY = 434,000 Ib/in.? A reasonable value for 
the yield strength of a plastic rigid material simulating mild steel would 
be Y = 70,000 lb/in.? giving a value of 6-2 for K, compared with 6-08 
from the theoretical analysis. 
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TABLE 1 





D (in. t (i.) 8D (in.*) P (lb) 





1°43 03 "006216 2,600 
1°43 *O4 ‘009572 4,260 
1°43 ‘ ‘013375 6,650 


"007569 3450 
“011652 4,99° 
‘016290 6,980 


‘OOS 451 3,700 
"013052 5,010 
"015250 7,450 
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; 
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3 6 wr 
Fig. 2. Plot of mean collapse load P vs. t!'VD from tests on mild steel tubes of 
various diameters and thicknesses, 


4. Discussion 

In view of the approximate nature of the analysis presented, it is 
surprising that there is such good agreement with experiment. The form 
of solution obtained agrees well, and the numerical magnitude is of the 
right order. The simple theory developed here does not take into account 
the effect of the superimposed axial stresses on the yield criterion, or 


detailed consideration of the deformation mode and equilibrium condi- 
tions. In effect, the methods used are those of upper bound techniques 
in limit analysis, which are not generally applicable to problems in which 
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large deformation or buckling occurs. No attempt has been made to 
determine the values of maxima and minima of the load during formation 
of the convolutions, since it is thought that a more realistic model of the 
true deformation mode would be necessary to give any realistic estimate. 

The main conclusion to be drawn from this analysis is that P = Ct'vD, 
where C is a constant best determined by a few experiments. 
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ON THE BUCKLING OF AN ANNULAR PLATE 
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SUMMARY 
This paper considers the buckling of an infinite plate supported along two con- 
centric circles and subjected to a uniform radial compression, or tension, along the 
inner circle. The solution is also applicable to a similarly loaded finite annular plate 
if there is a member of the requisite tensile stiffness supporting the outer circle. 
The effect of regularly spaced diametral supporting members is also investigated. 


1. List of symbols 
E Young’s modulus. 
v Poisson's ratio, taken as 0-3 for numerical purposes. 
t plate thickness. 
D Et /{12(1—v*)}. 
polar coordinates. 
inner and outer radii. 
radial and hoop stresses (positive if compressive). 
radial stress at r = r, causing buckling. 
deflexion of plate normal to surface. 
buckling stress parameter, to,7r?/D. 
r/o. 
11/To- 
number of diametral nodal lines. 
number of diametral supporting members. 
radial variation of w (see equation (4)). 
constants. 
indices defined by equation (7) 
(B—1)!. 


x, A, % indices defined by equations (21), (22), (26). 


2. Introduction 

When a uniform radial compressive stress is applied to the cireum- 
ference of a circular hole in an infinite plate, the resulting radial stresses 
in the plate decay inversely as the square of the distance from the centre 
of the circle; there are also tensile hoop stresses of the same magnitude 
and varying in the same way. An identical stress distribution exists in 
a similarly loaded annular plate, of thickness t, say, and bounded by 


circles of radii ry and r,, if the outer circle is supported by a member 


{ Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960] 
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of section area r,t/(1+yv). Such an annular plate will buckle when the 
applied radial stress reaches a certain critical value, depending on the 
boundary conditions. If the radial stress is compressive the mode of 
buckling will have rotational symmetry, but if the radial stress is tensile 
the plate will buckle into a number of circumferential waves due to the 
fact that the hoop stresses are now compressive. The present paper 
determines the magnitude of such critical compressive and tensile radial 
stresses and also considers the influence of regularly spaced diametral 
supporting members. The analysis throws light on certain buckling pheno- 
mena due to thermal or manufacturing stresses and on the buckling 
behaviour of highly tapered panels. 

The buckling of an annulus compressed along one boundary and free 
from stress along the other has been considered by Meissner (1), and the 
buckling of a uniformly compressed annulus by Yamaki (2). 


3. Differential equation of the buckled plate 
The differential equation of the buckled plate may be written in the 


form 


Weete teed a at) = ~ (1) 
r \r or r* of 


where compressive stresses are regarded as positive. Now the radial and 
hoop stresses in the plate are given by 


9 
ae ro\~ 
o, = —09 = |=") a9 


so that we finally obtain , 
i 1é 1 @\(@w lew 1 &w\ Blew 
poe ee ‘ea ot le re 
Or?  rér. r? 064]\ or® ) 


where the buckling stress parameter 8 is given by 


B = toyr?/D. 
The solution of equation (3) can be expressed in the form 
w = f(r)cosné (4) 
where n represents the number of diametral nodal lines and f(r) satisfies 
the equation 
d*f | 2 df _ (i+ 2n?2- —B) df. (1+-2n?- (1+ 2n?—B) yf a se . | 0 
di *rde Pr dF’ r’ ins 
(5) 


The particular case of buckling with rotational symmetry can be ob- 
tained by taking n equal to zero in equation (5). 
5092.49 C 
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3.1. Integration of the differential equation 
Equation (5) is homogeneous in r and the solution may be expressed in 


the form 4 
f= > Ap (6) 
i=1 
where p=T/To 
and, after some simplification, 
y? = 1+n?—4}8+(4n?—2Bn?-+ }f?)!. (7) 


When there is rotational symmetry n is zero and, as f is generally greater 
than unity, the indices y; can be expressed in the simple form 


Vi — +1, tid 
where ¢ = (B—1)# 


3.2. Boundary conditions 


(8) 


The following different conditions along the circular boundaries are 
considered: 


(i) zero displacement, 


(iil 


(ii) zero radial slope, 
) 


zero radial moment, 

(iv) zero shear resultant (i.e. the shear in the plate equals the component 
of the applied load normal to the deflected plate). 

Expressed in terms of f(r), these conditions are 

(i) f= 9, 

7° 
( 2f Vv d F L » 
aatea ant = 
df 1d*f 1 df  n*(3—v) 


4 ee a oe f  n(s—v) po. 
> aa” ak ie ali PY ar r pom 


(ii) 


(ili) 


(iv) 





4. Buckling with rotational symmetry 
We consider here the simplest case of buckling, with rotational sym- 
metry, which occurs when gy is positive (i.e. compressive) and there are 
no diametral supporting members. From equations (6) and (8) we can 
express f (or w) in the form: 
f = A,+A,p?+A,psin(¢ ln p)+A,p cos(¢ In p). (10) 
For any particular buckling problem there are four boundary conditions 
from which four equations connecting the relative magnitudes of A,, A,, 
Ay, and A, may be obtained; an expression for the buckling stress is found 
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from the vanishing of the determinantal equation. For example, if the 


inner and outer circular boundaries are clamped, conditions (i) and (ii) 
apply ; whence, writing p = 7,/To, 


(u2—1)(¢?— 1)sin(¢ In p) — 246{2p—(w? + 1)cos(Plny)} = 0. (11) 


Fic. 1 


Similarly, if the boundaries are simply supported, conditions (i) and (iii) 
apply, whence 


P - ¢?+ 1\?)_. Od fs ha 
(u?—1)i¢*—1-+ ( - pent nn) ~N—or + 1)cos(¢1n u)} = 0. 


| l+y 


(12) 
If the boundaries are free from rotation but are free to move axially, 
conditions (ii) and (iv) apply, whence 


p= 14 (7. (13) 
The buckling mode for this ‘freely clamped’ case is given by 

perforata] os 

If the boundaries are completely free, conditions (iii) and (iv) apply, 

whence 2 an Gg (15) 

and the buckling stress is thus independent of r,. This surprising result 

can be explained by a consideration of the buckling mode which is given 


by fo p—1. (16) 
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Substituting into equation (9, iii) it will be seen that the radial bending 
moment is zero everywhere; furthermore, from simple equilibrium, if there 
is no shear resultant at either boundary the shear resultant is zero every- 
where. It follows that the precise location of the outer boundary is, for 
this ‘free-free’ case, immaterial. 


Compressive buckling of annulus 


Fic. 2 


If conditions (ii) and (iv) apply along the inner boundary and conditions 
(iii) and (iv) along the outer boundary, we find 


¢ cos(d In x)-+-vsin(dln pp) = 0, (17) 


while if (iii) and (iv) apply to the inner boundary and (ii) and (iv) to the 
outer 


d cos(d In p) —v sin(d In p) 0, (18) 

Phe variation of the buckling stress with » for the six cases considered 

in this section is plotted in Fig. 2. The ordinate in Fig. 2 has been taken 
to be B(z—1)* uw, rather than £, for ease of presentation. 


4.1. Buckling stresses for infinite plate 
As r, tends to infinity the buckling stress depends only on the inner 
boundary conditions (ii) or (iii); if there is zero radial slope 


I t\2 
ee 19 
0” 12(1— v2) 4 (19) 
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and if there is zero radial bending moment 


EB (*. (20) 
12\r, 


5. Buckling with diametral nodal lines 

If there are m regularly spaced diametral members supporting the plate 
and a, is compressive, the buckling mode is of the form given by equation 
(4) with n equal to m. If o, is tensile the buckling mode exhibits diametral 
nodal lines naturally, i.e. even if there are no supporting members. 


saateeartnnggincaiaanianiaens es . —-y en 
| m=6 | 


, Circular boundaries clamped 





4 











Circular boundaries simply-supported 


a ae = 
2 3 a 5 


Compressive buckling of annulus supported on 6 diameters 
Fic. 3 


5.1. Compressive radial stresses 

Consider first the case of o, compressive (i.e. positive). For this case 
it can be shown that the indices y,; given by equation (7) are all imaginary 
and accordingly it is convenient to introduce the following symbols which 
are real and positive: 

x = {48 —1—n?+- (4n?— 2Bn? + }8?)4}4 (21) 

and A = {4B—1—n*—(4n?—2Bn? + }f2)h 14, 22 
rr . . AJ 
The buckling stress may now be determined by standard methods. For 
example, if the inner and outer circular boundaries are clamped the 
buckling stress is determined by 


(x?+-A®)sin(x In .)sin(A In x) —2 yA{ 1 —cos(y In p)cos(A In p)} = 0. 23) 
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Similarly, if the boundaries are simply supported 
1 : . - 2 2 ’ : 
x? +-A*- (* jsin( In .)sin(A In y.) - 
Ly 
~2yA{1—cos(x In »)cos(A In p)} 0. (24) 

The variation of the buckling stress with » for these two cases is plotted 
in Fig. 3 for the case when n (= m) = 6. 

As p tends to infinity it can be shown that 8 is determined by the 
vanishing of the term in parentheses in equations (21) and (22), whence, 
for m ia, Em? t 2 | 4) " 

09> - = 1+4-( 1 ——}°*}. (25) 
3(1—v?)\ro/ | m?} | 

5.2. Tensile radial stresses 

When a, is tensile (i.e. negative), it can be shown that two of the indices 

; given by equation (7) are imaginary and two are real. Because of this 


ee 


a srpreemnedseenenprousnements 
1) = ] te te — 5 a ——_— 4. ——_-—_ ++ 


5 7 8 
M 
Tensile buckling of annulus 
Fic. 4 


we introduce the symbols y and ys which are real and positive, y being 
given by equation (21) and & by 
yb = {1+n?—4}B+ (4n?—2Bn?+ 1B?) (26) 
If the inner and outer boundaries are clamped the buckling stress is 
then given by the equation 
(x?—-f?)sin(y In «)sinh(y In p)—2yf1 —cos(x In p)cosh (ys Iny)} = 0. 
(27) 
If the inner and outer boundaries are simply-supported the buckling stress 
is given by the equation 


721 ys? 2 : ' 
pew ww | sin( In y)sinh (ys In p)— 


—2y%{1—cos(x In p)cosh(yInu)} = 0. (28) 
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If there are no diametral supporting members the method of solution 
of equations (27) and (28) is to assume various integral values of n for 
each value of y, and solve for 8; the envelope of the resulting family of 
curves gives the true value of 8 (and n). Fig. 4 shows the result of such 
calculations. As pu tends to infinity it can be shown that for both clamped 
and simply supported cases 


si —E [t) 
™ °0* 4(1—)\rq 


A comparison of this result with the corresponding result for compressive 
do, &.g. equation (19), shows that the buckling tensile stress is three times 
the buckling compressive stress. 

If there are m regularly spaced diametral supporting members the 
method of solution of equations (27) and (28) is to assume, for each value 
of 4, various values of n which are integral multiples of m. The envelope 
of the resulting family of curves gives the true value of 8 (and n). As u 
tends to infinity it can be shown that, for m > 2, 


sia oi iat). (30) 


~ 12(1—2) \r, 


(29) 


6. Conclusions 

This paper considers the buckling of an annular plate in which the radial 
and hoop stresses vary inversely as the square of the distance from the 
centre. The presence of diametral nodal lines, either imposed or occurring 
naturally, has been allowed for, and this has enabled the buckling 
behaviour under tensile, as well as compressive, radial stresses to be 
determined. Figures are presented giving the buckling stress as a func- 
tion of the geometrical proportions of the plate for a variety of boundary 
conditions. 

The class of problems considered is unusual in that exact solutions are 
available in terms of elementary functions despite the fact that the stresses 
are varying throughout the plate. 
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THE THIN-WALLED CIRCULAR CYLINDER SUBJECTED 
TO CONCENTRATED RADIAL LOADS 


By L. 8. D. MORLEY 
(Royal Aircraft Establishment, Farnborought) 
| Received 7 May 1959] 


SUMMARY 
Using a new equation governing the radial displacement of a thin-walled circular 
cylinder, expressions are obtained for this displacement and the resultants in an 
infinitely long cylinder subjected to equal and opposite concentrated radial loads. 
In addition, a method of numerical approximation is presented which enables the 
calculations to be carried out with rapidity and yet retaining an accuracy usually 
associated with Fligge’s equations, 


1. Introduction 
ProBLeMs of equilibrium of thin-walled circular cylinders under various 
loading conditions have attracted the attention of many authors. The 
solution of some of these problems has required quite complex analyses, 
but it is significant that when numerical values are needed, resort is usually 
made, notwithstanding loss in accuracy, to simple approximations so that 
the computations are not prohibitively laborious. 

For example, during the last few years, Hoff (1, 2) and his collaborators 
have employed a simple approximation which was introduced by Donnell 
(3) in 1933 whilst investigating the torsional stability of thin-walled 


cylinders. Unfortunately, however, the error in this approximation 
becomes intolerable beyond a certain wavelength of circumferential distor- 
tion and, recognizing this, Hoff (4) later examined the accuracy of the 
trigonometric eigenfunctions arising from the equation which governs the 


behaviour of the radial displacement. For this purpose a more accurate, 
and more complex, equation derived by Fliigge (5) was used as the standard 
for comparison and Hoff found that it was necessary to restrict the range 
of application of his previous equilibrium solutions. 

Earlier, in 1946, Yuan (6) had used the same Donnell equation to 
determine the radial displacements in an infinitely long and thin-walled 
circular cylinder subjected to equal and opposite concentrated radial loads. 
However, this equation is particularly inaccurate for this purpose and so, 
in 1957, Yuan and Ting (7) presented another solution to the same problem 
but this time using the more accurate Fliigge equation. This second 


+ Paper prepared whilst attending the ‘ Post Graduate Course in Structures and Materials’, 
Cambridge University. 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960] 





THE THIN-WALLED CIRCULAR CYLINDER 25 


solution is much more complex but, although he found the computation 
tedious, Yuan discovered that his first solution gave a result some 25 per 
cent in error for a particular numerical example. 

More recently a new equation has been introduced (8) and, although 
it has the simplicity associated with Donnell’s equation, it has been 
demonstrated that its trigonometric eigenfunctions are in close agreement 
with those calculated from Fliigge’s equation. Now the purpose of this 
paper is to present an application of the new equation to the solution of 
Yuan's (6, 7) equilibrium problem which has just been described. The 
analysis follows that employed by Yuan but, in addition to presenting 
an expression for the radial displacement, expressions are here given for 
the stress and moment resultants. It is further shown that important 
simplifications can be made to these expressions which enable the numerical 
calculations to be carried out with rapidity and accuracy. In the example 
it is shown that a brief slide rule calculation yields a value for the maxi- 
mum radial displacement which is only 0-5 per cent different from that 
obtained by Yuan (7) in his second and accurate solution. The method 
of carrying out these simplifications is not confined to the present problem, 
but provides a basis for the ready and accurate solution of many further 
problems of thin-walled circular cylinders. 

Finally, Fligge’s equations are used for the stress and moment resultants 


although, as Hildebrand, Reissner, and Thomas (9) have pointed out, they 


contain terms of small order whose retention is not justified by the basic 
assumptions. However, the inclusion of such terms does not reduce the 
accuracy, and it is a simple matter to discard them as and when the 
numerical analysis shows that this is justified. 


2. Notation 
a mean radius of the thin-walled cylinder. 
D flexural rigidity, D = Eh3/12(1—v?). 
E Young’s modulus of elasticity. 
h wall thickness. 
i AT 
K non-dimensional structural parameter, 4K4 = 12(1—v*)(a/h)?. 
m first even integer satisfying 5m?(2m?—1)* > 16K4. 
M, bending moment resultant along a generator. 
My bending moment resultant in the circumferential direction. 
M9, Mp, twisting moment resultants. 
n even integer. 
direct stress resultant along a generator. 
direct stress resultant in the circumferential direction. 
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shear stress resultants. 
concentrated radial load applied at 2 = 0; 6 = 0, 7. 
distributed radial load. 
shear stress resultants acting normal to the surface. 
displacements respectively along a generator, a circumference, 
and a radius. 
non-dimensional displacements defined by u = u'/a, v = 
w w/a. 
distance measured along a generator. 
non-dimensional distance defined by x = z’/a. 
(2n?—1)/4K2. 
1+ K4 nt. 
introduced by equation (1). 
roots of the auxiliary equation 
(A?-+-n?)?(A?+-n?— 1)? +- 4K 44 = 0. 
Poisson’s ratio. 
angular coordinate. 
Subscripts following a comma indicate differentiation, e.g. 


Ww 299 = Cw/éx?06?. 


Laplace's operator is denoted by V? (= é?/éx?+-é?/0?). 


3. The radial displacement 

The thin-walled ciretlar cylinder subjected to equal and opposite 
radial loads is shown in Fig. | and the positive directions of the displace- 
ments and resultants are given in this figure and Fig. 2. 

The method of analysis follows that employed by Yuan (6) who began 
by replacing the concentrated radial load P by a distributed load q ex- 
pressed as a function of the coordinates along a generator and a circum- 
ference. Thus, 

: jsin(Ad a) 9 S sin(nd a)sin(Ad /a) mi = 


(1) 
Now, unlike Yuan who used first Donnell’s and then Fliigge’s equation to 
govern the behaviour of the non-dimensional radial displacement w, the 
new equation (8) 


FP .. 
q(x. @) = lim oe oho esate: ul I cen 
aa ma? 5.9 J | (Ada) pe (nd/a)(Ad/a) | 
s n=2,4,... 


a> 
a 
5 VG 


IIL I 


V4(V2+-1)?w+4K tw 


is used. Substituting equation (1) into (2) and taking the limit as 6 
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approaches zero yields 


2a Pi l bf aoe cos(Az) 


teenie m®D\2 J (®—1)?+4K4 


dxr+- 





a r (A? +n?) cos(Az) 
p3 cos n6 | OR Ln Ea . (3) 
n= <,4,... 0 


Fic. 2 


The integrals in this equation are readily evaluated with the aid of the 
calculus of residues and it is found that 


w(x, 0) = Sl Pca x(K®—})¥}{(K® + })teos[a( K?-+$)#) + 


+ (K?—})sin[a( K?+ 4)*}}— 





—2 S in| Ain (Agn +? )etrrn Aon(Adn + n2)ettrrn 
~ (Ain +n? — 1)(Ajn— n+?) © (AR, +0? — 1)(03,—n4 +n?) 
1 2n 


cos no 
(4) 


n=2,4,. 
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where the A,,, and A,,, are the roots of the auxiliary equation 


2n 


(A2-+-?)?(A?-+-n2— 1)? + 4K4A4 = 0. (5) 


Referring to the previous paper (8), the roots of equation (5) are given by 


2 211 3K i(2n2—1))4 
Ain —n*-+-} ike) |) 


oe | (6) 
n2-+}4 iK?|1 ate | 


| ro | 
- 


The above equations are noticeably simpler than the corresponding ones 
derived by Yuan (7) in his second and more accurate analysis employing 
Fliigge’s equation. Nevertheless, a close inspection reveals the possibility 
of further simplification which enables the numerical calculations to be 
carried out with satisfying rapidity and accuracy. 


4. A method of numerical approximation 

\s noted by Hoff (4), the non-dimensional structural parameter A must 
have a numerical value greater than 5 because for lower values the assump- 
tions of thin shell theory are invalid. On the other hand, A must be less 
than 50 otherwise the shell is so thin-walled that it can hardly fulfil 
structural requirements. 

Bearing these values in mind and making the substitution 

2n?— 1 
4k? 
then the first of equations (6) can be rewritten as 
Mf, = K-28, —i{1—(1—4i8,)}} 
(1 215, + 28? + 4183 —1084 t .-)}] 


21K *62(1+-275,,), when 582 is neglected in comparison with unity. 


In fact, the roots A,,, and A,,, of equation (6) can be approximated by 


) a Ay, = (1+i)K(1+28,), (8) 


provided that it is agreed to neglect 552 in comparison with unity.f 


BY 


+ With regard to the accuracy, it is interesting to note the following. Fligge’s auxiliary 


equation 1s ) 2\2 5 2 ‘a 
! (Ap-+n*)2(AZ, + n?— 1)? + 4K4A$4 2(1—v)(AZ—nt+n?*) = 0, 


so that if A is a root of the present auxiliary equation (5) then 
Ap (1+ e)A, 

where the error ¢ is approximately 
—(1l—v)(A*—n*+n?) 


2 2(? +n? — 9) (0? +n? — 1)? +n?) + 4K)? + (1—v)(3A*— 8 + n?)] 


¢€ a 





from Newton's approximation. If now the rough approximation is made of neglecting 43, 
in comparison with unity, it is found on employing equations (8) to (10) that the error ¢ 
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Furthermore, within the same approximation, 
Ain—n*-+n? = —n*(n?—1)) 


Ai, —n*+n? = Af, J 


(9) 


The accuracy of equations (8) and (9) is at its best for the terms containing 
the smaller values of n and these terms are usually the most important 
constituents of the final solution. 

Employing the approximations (8) and (9) together with the identities 


2 5 9 - ‘ 25 Rr 2\2 
(Az, + n*)(Aj,, +n? — 1) —2iK*Aj,, ) 


N (10) 
(AS, +n®)(A3, +-n?—1) = 2iK203, J 


which are valid for all values of n, the relevant terms of equation (4) 
simplify to 


Ay (AZ, +2?) (l+t) ( n? | (Poe)! 
4k2- 


— dy 


(AR, +n? 1)(Ah,— attr?) K(2n?—1)|n®—1 


Az, (AZ, +n?) . (+i) 


| (11) 


(A3,, +-n?—1)(Ad, —nt+n?) ~ 4hk% 

The above approximation is not sufficiently accurate for the larger 
values of the even integer n, but it is then permissible to neglect unity in 
comparison with n*. If, furthermore, it is assumed that 

A?+-n?—1 == A®+-n? 
then the roots A,,, and A,,, of equation (6) can be approximated by 
Ain = K—nv2(«,—1)#+ t{nv2(«, + 1)#—FK} | 
2r,,, = K+nv2(x, —1)'+ifnv2(x,+1)!+ K}) 
' Kt 
where ; 1+—, 
4n* 


is nearly independent of n and 

l—y 

4hk?° 
Therefore, consistently with this value of ¢, it is sometimes permissible to neglect quantities 
greater than 582 in comparison with unity. On such occasions, C(n, K)32 can be neglected 
in comparison with unity, where 


€ 


. 1 /l—p’ 
Cin, K) ws (axe)- 
Some typical values of C(n, K) are tabulated below for Poisson’s ratio vy = 0-3: 
n y 2 4 2 4 6 
K t 10 10 50 50 50 
C(n, K) <i 5-7 5 143 73 <5 


Finally, since 5, > 10 there is justification in retaining the term 3, in the expression 
Aen = (14+ i)K(14-48,). 
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and the identities of equation (10) can be replaced by 
Az, tn® = (1—i)KAj, ) 


14 
dB +n? = (1+i)Kp, | 6) 


Equations (12) to (14) are identical with those obtained directly from 
Donnell’s equation and their accuracy improves as the value of n increases. 

Using equations (12) to (14) the approximation corresponding to equa- 
tion (11) above is 


= —_ AinAtn +n?) ilies i 

(A? +n?—1)(A$,—n*+n?) © 2K{KA,,—(1+i)n?} 
sii Aan (Abn +n#) _ — | 
(A3,,-+n?—1)(A$, —n*+n?) © 2K{KA,,—(1—i)n?} 





(15) 








where the roots A,,, and A,,, are obtained from equation (12). 
In a numerical application this second approximation is superior to the 
first one given in equations (8), (9), and (11) whenever 


5n*%52 > 1, 
or, on substituting from equation (7), whenever 
5n?(2n2?—1)? > 16K4. (16) 
In the special case when z is zero, equation (4) for the non-dimensional 


radial displacement simplifies to 


I 





(F) (sr). : K _ 2K 7 2K*n? 


| cos nO+- 


h an om Ly | (n?—1)(2n?—1) 


~ : 
(nt 1) co8n8, (17) 
nx, 

n= m,m + 2,... 


where m is the first even integer satisfying 
5m?(2m?—1)? > 16K*. (18) 


The second summation in equation (17) was first derived by Yuan (6). 

The errors in the physical quantities, such as the displacements and the 
resultants, are of the order (h/a)x 100 per cent and, consistently with 
this, it is assumed in the sequel that 


ha is negligible compared with unity (19) 


for the first term and for those containing the smaller values of n. This 
has already been done in equation (17). 
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5. Numerical results 


The maximum radial displacement occurs underneath the point of 
application of the concentrated radial load and Yuan (7) has calculated 
the values of the displacement quantity (a/h)(Zh?/P)w(0,0) using both 
Fliigge’s and Donnell’s equations. His results are reproduced in Table 1 
below along with those obtained from the present equation (17). The 
calculations are appropriate to a thin-walled cylinder with (a/h) = 100 and 
a Poisson’s ratio v = 0-3, which gives a value for the non-dimensional 
structural parameter of K = 12-85. The value of the even integer m is 
soon found from equation (18) to be m = 6. 


TABLE 1 


Comparison of the Fourier coefficients for the displacement quantity 
(a/h)( Eh?/P)w(0, 0) 





n 2 4 6 8 10 





Present equation (17) 52 101 43 24 14 
Fliigge . : ; - 5: 99 44 24 14 
Donnell . : , 348 91 43 24 14 























It is seen from Table 1 that there is little difference between the results 
obtained from the present simple calculation and those obtained at length 
by Yuan using the more accurate but very complex Fliigge equation. 
Furthermore, whereas Donnell’s equation yields a value for the displace- 
ment w, underneath the load, which is 25 per cent in error, the present 
accuracy is to within 0-5 per cent.t 


6. The resultants in terms of the displacements 

Attention is confined to the bending moment resultants My, M, and the 
direct stress resultants Ng, N, because it seems that these are the quantities 
most commonly required in practice. Lf needed, however, the twisting 
moments M,», Mg, and the shear stress resultants Ni», Np,, Q,, and Q» 
can be dealt with in the same way. 

It is customary and concise to express the resultants in terms of the three 
displacements u, v, and w, but, for their calculation, it is usually more 
convenient for them to be expressed in terms only of the radial displace- 
ment w. 

Now, when u and v are eliminated from Fliigge’s resultants (which are 
listed in the previous paper (8), they become 


) 
M, — _- (w gg-+w-+v0 ,,}, (20) 
tf 


t See Appendix. 
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1 fh\*( 0) op ywew ) ~3(l—v) 1/h\2 7, 
a al) ‘a 2v)V2w 700 2 ti : ._ wats) W xxez00) (21) 
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hoc tv(w gg+-w)}— (1 —v?)w 2269 + 


Ving = a LV4(V2-+ lw 
a2| ere: 


—v)( Ww 9909 — a aie 


3x(1—v) 1 ( 


2 12\a 


F | 
] W rrrx00 | 


-aVi4q i Pieove. tye (1—v)( 
ri 


: : — 
| 0000 — © zrre T =U 00),22+ 


3v(1—v) 1 /h\* ) 
i pea 12 ‘) eed 


) ¢ , 2 
VAN, 4 K4w, roo — 2(1 v)V?w 206 —— " (-) Wzazstt) (23) 





a*| y l2\a 


These equations contain many terms which are so small that their reten- 
tion is not justified by the initial assumptions. When the substitution is 
made for w from equation (4), and the numerical approximations described 
in section 4 are applied, these small terms are revealed and can then be 
discarded. However, before writing down these final equations it is 
interesting to make a premature application of the approximations and 
so rewrite the above equations retaining only the significant terms. They 


are D 


M, = (W gg tutu, (24) 


a 


D, 


{uw .,+v(wagt+w)}, 
a . . 


or N, —aq + Pouye L1)w, 
a* 
4DK4 


5 “ rr00- 
a* 


In arriving at equations (24) to (27)t it is necessary to neglect some terms 
+ The corresponding equations obtained from Donnell’s approximation are 
D 
Mo — — (weet vw xz), M, “2 (Ww + VW py), 
4DK* 
’ w 


a- 


: D 
ai 2x72 or Ne —aq-+ aa ts 
4DK* 


a? 


vi N, WU, -200* 
When q and w are independent of x then the equations should satisfy those defining the 
behaviour of a circular ring, i.e 
l : i dad .\3 
No+ = Me.e3 Nope + Ne = —29 and (ss + 1) Ww = 


Equations (24) to (27) satisfy these, but Donnell’s equations do not. 
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which are only noticeable when the wavelength of distortion is so short 
that the assumptions of thin shell theory are invalid. The second of 
equations (26) reveals the expected direct relationship between circum- 
ferential stress resultant and applied radial loading. 


7. Calculation of the stress resultants 
Starting with the circumferential bending moment M, and substituting 


for w from equation (4) into equation (24) there results 
P [ve-** . . 
M,(x,0) = —= oR (sin xK —cosxrk)- 


“TT 


{A,,( (Ain +n?—1)(AZ, +n? _ 


—— -+- 





2 +-n?—1)(At, —n*+-n?) 


LA (vAz,, +-n?—] )(AZ,, + n2)e'tAan) 


en ‘ 
(Az, TNA At —nt in) joo nd). (28) 
where equation (19) is used to simplify the first term. Now, before apply- 
ing the numerical approximation described in section 4, it is advantageous 
to take note of the consequences of equation (19) for the smaller values 
of the even integer n. The consequences are, for Poisson’s ratio v = 0-3, 
that 1 

P 3K: 8 negligible compared with 0-3 


(29) 


2 
In 


is negligible compared with vor | 


1 is negligible compared with 0-3 Aj, 
Using these, together with equations (7) to (10), the relevant terms in 
equation (28) simplify for the smaller values of n to 


Ain (vAin Lm? “Win +n?) | (l+i)n?  (1—i)(1-4 ) 


(A2, + n?—1)($,—nttn®) © K(2n?—1)' 4K 
Ay (vA$, +-m®—1)(AB, +n?) (Li) (14-4)(1+v)n? | 


: (Az, +n? Das — a in?) 2K 4Kk3 


(30) 


This approximation is insufficiently accurate for the larger values of n, 
but it is then possible to neglect unity in comparison with n®. Using 
equations (10) and (12) to (15), equation (30) becomes 
Aan(vAfntn?— 1), +m). (Li) i(1-+-v)n2 
(3, +-n2—1)(At,—nt+n?) 2K 2K{KA,,,—(1+1)n*} 
Az, (vAz, +n? —1)(A3, +2). v(1—i) i(1-+-v)n2 








, (31) 


(Az, +-n2—1)(A$ i. —nitn?) ~ 2K  2Ki (KA, —(1 i)n*} 
where the roots A,,, and A,,, are given by equation (12). Equations (28), 
(30), and (31) combine to give a particularly simple expression for My when 


5002.49 D 
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x is zero, i.e. 


a Je * | 2n? vy n(1+y)) 
M08) = slot D2 \K@ay tet me jenet 


ae 2 44) 
4 VA +) > (Kn +1) cos nd], (32) 


» 
“ nk, 
n= ™M,M+ 2,.0. 


where «,, is given by equation (13) and the value of the even integer m is 
inferred from equation (18). The infinite series converges as cos nin 
which is too slow for practical computation, but advantage can be taken: 
of the known summation 


SY cosnl _ _ 4 log(2sin 8). (33) 
wien. 

Underneath the point of application of the concentrated radial load, where 
both x and @ are zero, it is seen from equation (32) that the bending moment 
is infinite and this is in agreement with the results of the elementary theory 
of bending of flat plates. 

The explicit equations are not given for the bending moment resultant 
M, along a generator because, from equations (24) and (25), it is readily 


seen that (Me My, a | 
vel 


34 
D D (34) 


= (1+) 
The stress resultant Vy in the circumferential direction is given by 
equation (26) and on substitution for w from equation (4) it is found that 


> 
aq— : Ke *K(cos xK +-sinak)+- 
7a 


J 9 9\4 9 ’ 9\9 = 
> im! 42 n(Ajn +n?) path Aan (Aan +n?) ia loos né (35) 
—_ (At, —n*-+n?) (Ai,,—nt+n?) | 
N= 2,3,.06 " 


The relevant terms in this equation can be approximated by, for the 
smaller values of n. 


Ai, (AZ, +2)? . = (1+)n? 


A,,, (Az, +-n?)? (1— a | 


(36) 


Sanl@en TR! oe (144)K 4+ 
(A3,—n*+n?) © devi 


2K 


where equations (7) to (9) and (29) have been used. For the larger values 
of n, where it is possible to neglect unity in comparison with n?, equations 
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(36) become 
Ay (Ain- +n" ‘asta 
n\ln = K + 
(At a —n*+n?) * +o 


Aen (An + 0?)® = (1+i)K + 


(Aj, —n*+-n?) 


nk 
' KM — a Tine | 
a ) 
KA,,—(1—i)n? 

2). When 2 is zero, 


am). 


(37) 


where the roots A,,, and A,,, are given by equation (1 
equations (35) to (37) combine to give 


P m—2 
N,(0, 0) = —aq— ‘ (K?—n?)cos né + 

2ra\ 
n=2,4,... 


oa 
+K*y\2 = cos nd, (38) 
nk, J 


n=m,m-+ 2,... 


where «,, is given by equation (13) and the value of the even integer m is 
inferred from equation (18). The infinite series converges quite rapidly. 
Finally, the stress resultant \, along a generator is given by equation 


(27) which, on substitution for w from equation (4), becomes 
] 


> D 7423 pirrin 
NAz, 0) = : S im/!- - m2. Min : ae" te 
ma \(A2, +n? 2— 1)(A2, + n?)(At, —nt-+ n?) 


4K4n?)3, eithen | : 
+ OE FETE, + nF, nt pw] OME OO) 





For the smaller values of n, the relevant terms in this equation can be 


simplified to 
_ —— ss ; _ K(2n* - 1)(1 1 a) 
OF aE Fat atta) = — Bey — |! + aR) | 
‘4K n2X3, - _ n(1—i) | 
—nttmn?) ° 2K 


(Az, + n2)(A3, +n?—1)(As, 


(40) 


where equations (7) to (10) and (29) have been used. For large values of 


n, where it is possible to neglect unity in comparison with n*, equation 
(40) becomes 
4K nr}, 
(A2. 4 n2—I1)(at ww Ge i 
n*)(AZ,, +n?—1)(At,, —n4 +n?) Ain- 


Ain 
4K4n?A}, 





eK 
-(14 + -i)n? | 
n?kK 


OEM )OR Fm AR — wpa) ~~ KA, — ant 


(41) 


where the roots A,, and A,,, are given in equation (12). When z is zero, 
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equations (39) to (41) combine to give 


> m—2 oe 2 
N.(0, @) = P a \(<— ge OC cs n+ 
P 2ra| Ky \\n?—1 2 | 


x 
a K,—1)# 
+ K2y2 > («n—}) cosnO}, (42) 
nk, 
n=m,mM + 2,... 


where «,, is given by equation (13) and the value of the even integer m is 
inferred from equation (18). The infinite series in this last equation is 
identical to that in equation (38) for the stress resultant \,(0,@) in the 
circumferential direction. 

No numerical values are given for the resultants because it seems that 
there are no values determined experimentally or calculated from Fliigge’s 
equations with which to obtain a comparison. It has, however, been 
demonstrated how to obtain expressions for the resultants which are 
simple and which have an accuracy of the order of (/a) x 100 per cent. 


APPENDIX 
In a further and recent paper, which came to the author’s notice after completion 
of the present paper, Ting and Yuan (10) have again re-examined the problem but 
this time using the ‘complete Donnell equation’ 
a? 
vs edhe V‘q, 


. Das " "Says a 
V8w + 2W ggegeg + W og9g + 4KAW pre = D 


which was introduced by Donnell in 1938 (11). This equation, which contains all 
the vital constituents present in Flugge’s equation, seems to have remained un- 
known and unused until now. 

Corresponding to equation (17), Ting and Yuan obtain 


h? 2K 4 K4v2 — (Kk, +1)4 


é ) 
rents, OS nO — / a n, 
n-2)i 7 er nk, 
n five n=m,m+ 2,... 


= }w(0, A) 


77 


where m is now, presumably, the first even integer satisfying m? > K. 
For ash 100, this leads to a value of m 4 and the following supplement to 


Tabl ] above. 
SUPPLEMENT TO TABLE 1 


Fourier coefficients for the displacement quantity 


(a/h){Bh?/P)w(0, 0) 





n 0 | 2 





‘complete Donnell | 
equation’ 520 91 
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NOTE ON A PAPER BY G. M. L. GLADWELL 
‘SOME MIXED BOUNDARY VALUE PROBLEMS OF 
AEOLOTROPIC THIN PLATE THEORY’ 


By F. J. HAWLEY 
(City of Portsmouth College of Technology) 
{Received 4 June 1959} 
SUMMARY 
In section 5 of Gladwell’s paper (1) it is stated: “The author believes, but has been 


unable to prove, that «,, Kk, are positive in the general case.’ This note proves that 
K,, kK, are positive for the general case of digonal symmetry. 


1. Introduction 


Gladwell (1) in his discussion of the mixed boundary-value problem for 
the half-plane introduces a parameter « satisfying the equation 


(44 — BC)x?+-( B?+C?—2AA)xe+(AA—BC) = 0 


where 


B 





‘al 
and Yi; Yo are the two roots of 


a,+ 2a, y+ (2a,+a,)y?+ 24,7 +4,y' = 0 (1.2) 


satisfying y, <1, y. < 1. Gladwell proves that the roots x,, x, are 
always real and that they are positive when the material is orthotropic. 
It will now be proved that x,, «x, are positive in the general case of digonal 
symmetry; this completes Gladwell’s solution of the mixed boundary- 
value problem for the half-plane. 


2. Method 
If we put y, rei, y, = re!’ where r, < 1, r, < 1 and write 
mtr, 8 =—7f,+I1/r, 

then it may be shown by using the relationships between the roots and 
coefficients of (1.2) that 

Ret 6+), 

—} Re'(s,+ 8, e'?), 

2a, R(s, 8.+2cos¢) 


{ Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960) 
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so that 


—" l R? B.. By 
BC—AA =a, Rr, yt | — |a+rra(a+3)+ 


i--> 
i°2 4 


sseng-dimd)p en 


The positive-definiteness of the elastic energy function demands, in the 
general case of digonal symmetry, that the determinant 

a, a, ay; 

a, 

@y G, 4, 
and all its principal minors shall be positive. The conditions a, > 0, 
aj > a34,, and G > 0 yield 


9 
- 


} R21 +r} 7r3)(1/ri+ 1/73) < a, Rs, 8.—2a} 


and a, > Reos¢. 
If (2.2) is substituted into (2.1) it is found that 
BC—AA > }R(a,— Reos¢)(r7y—1/7,)(r72— 1) 
>0 


and this is the condition for «,, «, to be positive. 
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THE DAMPING EFFECT OF DISTRIBUTED AND 
CONCENTRATED RESISTANCES ON SMALL 
PERTURBATIONS IN A UNIFORM FLOW 


By A. P. BURGER?t and T. W. van DER LINGENT 
[Received 5 February 1959] 


SUMMARY 

The damping effect of a distributed resistance on small stationary harmonic 
perturbations in a uniform incompressible flow is theoretically analysed for the 
two-dimensional case, A quadratic friction law is assumed. 

A solution is obtained for the damping factor and simplified versions are found 
for the cases of widely distributed resistances and concentrated resistances. These 
solutions have applications to the damping effect of flat plate heaters used in blow- 
down tunnels and gauze wire screens respectively. It is shown that there is good 
agreement between the theoretical result for a concentrated resistance and experi- 
mental results achieved on gauze screens by previous investigators. 


Notation 
half-length of resistance plates. 
friction factor. 
drag coefficient. 
perturbation wavelength in y-direction. 
wave number = 27/1. 


perturbation velocity components. 


space coordinates. 

non-dimensional half-length of resistance plates. 

flow regions as indicated in Fig. 1. 

non-dimensional wave number. 

constant speed of basic current. 

total flow speed. 

non-dimensional space coordinates. 

integration constants. 

damping factor. 

vorticity. 

solutions of equation (6). 

(perturbation) pressure. 

fluid density. 

perturbation stream function. 

x-dependent part of w. 

non-dimensional form of F. 
+ National Physical Research Laboratory, South African C.S.1.R., Pretoria. 
t National Mechanical Engineering Research Institute, South African C.S.I.R., Pretoria. 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960] 
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1. Introduction 


The damping effect of concentrated resistances, such as gauze screens, 
on small perturbations in a uniform flow has been analysed by Prandtl (1), 
Collar (2), and Taylor and Batchelor (3). The results of the first two 
authors were embodied by Taylor and Batchelor in a relation for damping 
in which they used a factor « determined by the directions of the per- 
turbations when entering and leaving the screen. They established the 
value of « by means of experiments on an undisturbed flow and found good 
correlation with the experimental results of Collar and MacPhail (3) for 
gauze screens. 

Recent advances in the technique of pressure driven intermittent wind- 
tunnels have, however, brought for ward the need for an analytical method 
of determining the damping effect of a distributed resistance in a uniform 
stream. A flat plate temperature stabilizer of the type conceived by van 
Spiegel (4) may be used to damp out perturbations in the tunnel flow, if 
it is placed behind the blow-down valve. The temperature stabilizer 
consists of flat plates set parallel to each other and to the air stream 
which, by virtue of their heat capacity, limit the drop in temperature of 
the air stream during the blow-down period. 

With such an arrangement two possibilities, which may be treated two- 
dimensionally, arise. The first is that of perturbations distributed in a 
plane at right angles to the plates. If incompressible flow is assumed, the 
velocity between a pair of plates will be constant and the flow at the exit 
will be parallel to the plates. In (3) it is shown that in this case the 
original solution due to Prandtl will apply. 

The second possibility is that there are perturbations in a plane parallel 
to the plates, so that the flow between a pair of plates is to be considered 
and the plates themselves do not constrain the flow. This case is treated 
in the present analysis. If the static pressure at normal sections of the 
flow is then assumed constant, a useful approximation to the damping 
effect may be found in a simple manner. In section 5 it is derived in the 
form uu, = e-**¢, where u,/u, is the ratio of the amplitude of perturba- 
tions before and after damping and 2ha is the total resistance coefficient 
of the plates. 

From physical considerations it may be deduced that the condition 
of constant static pressure at normal sections will apply when the length 
of the element forming the distributed resistance is large compared with 
the scale of the perturbations. In this case sharp curvature of the streamlines 
with the attendant normal pressure gradients will not be required to 
redistribute the flow. 
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However, it is not clear from this result how long a resistance element 
should in fact be, before the above damping equation will apply. On the 
other hand, the solution for gauze screens is applicable when the length 


of resistance element is very short compared with the scale of the perturba- 


! 


tions that are being damped. 





L 


Vv. 


{| U+u, 











; 

‘ 
al * *) oO 

Fic. 1 

In order to determine the damping effect of distributed resistances, a 
solution must therefore be found in terms of the length of the resistance 
element and the wavelength of the perturbation. This solution may be 
expected to yield the results for a very long resistance element and for a 
gauze screen when the ratio of wavelength to length of resistance element 
tends to zero and to infinity respectively. The solution derived below 
vields the above-mentioned constant static pressure result, and a useful 
approximation to the behaviour of gauze screens, in the respective limiting 


cases. 


2. Idealized arrangement 

The idealized two-dimensional arrangement, in which a solution for the 
damping is sought, is shown in Fig. 1. We consider a uniform two- 
dimensional incompressible flow of velocity U in the X-direction, with 
small stationary perturbations u and v in the X- and Y-directions respec- 
tively. The perturbation passes over a region L from x = —« tox = —a 
in which there is no resistance, followed by a region M from x = —a to 
x +a in which resistance is proportional to the square of speed, and 
by a final region NV from « = +a to = +0 in which there is again no 
resistance. 

This corresponds to the set-up of Taylor and Batchelor (3), with the 
gauze screen replaced by a third flow region subject to distributed 
resistance, 
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3. Flow equations 

As the third dimension, normal to the resistance plates, has been 
omitted, the resistance is taken to act uniformly on the flow in the middle 
region M. We assume a quadratic law, so that the magnitude of the 
resistance force per unit mass is 


shV? 


where V is the speed and fh is a proportionality factor. Denoting the 
pressure by II, the equations of motion are then 


v= Me _ sayy’ (1) 
p 
fl 


Pp 


v_$hVo, 


in which dots denote time derivatives, subscripts x and y denote partial 
differentiation, while u and II, with primes indicate values pertaining to 
the total field of flow, i.e. basic current plus perturbation. 

Forming the curl of these equations, and denoting vorticity by ¢, yields 


{= Shi VE+eV,—(U+upV} 
which, on linearization with respect to the uniform basic current, becomes 
Ul, = —th(UC—UD,) = —jyhU(v,—u,—1,) 
i.e. C, = A(u,—4,). 
In the regions L and N of free flow. h vanishes, so that we have in the 


three regions 


a (h(u,— bv.) (a, <a) 


a (ja! >a). 


Since the flow is assumed incompressible, a perturbation stream func- 
tion & exists such that u = —y, and v = w,, and the equations become 


—_ (A(- byy— bez) (2| <a) 


Vip, = 


lo (jz| >a). 


If the problem is further simplified by making the variation of the 
perturbation sinusoidal in y, so that 


% = V(x)sin py 
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(hipey—4¥") (|x| <a) 


we get yr" — pt" = = 
lo (|x| >a) 


where accents denote derivatives with respect to x. Finally, non-dimen- 


sional quantities are introduced by taking //A as unit of length and U as 
unit of speed, viz. 


X,Y, A = hea, by, ha, 


i 


¢(X) = * (2), 


P 


h 
The equations then become 
}” +4¢"— P*(¢’+¢4) = 


$” — Pr’ 


4. Boundary conditions 

If the lateral walls are parallel to the basic current, their positions may 
be assumed to be situated at the zeros of sin py; the boundary condition 
v = 0 at the wall is then automatically satisfied. 

If it is assumed that the plates are infinitely thin, no concentrated 
forces act on the flow at their edges = -La, and we therefore require the 
stream-function and its derivatives of first and second order to be 
continuous, 

On the other hand, the blockage effect may be taken into account by 
following Prandtl’s discussion on vortices entering a contraction, as given 
by Uberoi (5). If the contracted section is designated by suffix c and 
the open section by suffix 0, and if the blockage factor is defined by 
U, = cl), the relations u, = (1/¢)u, v, = Vo, and f, = (1/c)f_ are obtained 
for the two-dimensional case, i.e. ¢, = cdo, $, = bo, $. = Cho. If these 
conditions are used instead of the above-mentioned continuity conditions, 
the damping factor obtained is not influenced, to first-order approxima- 
tion, in the two limiting cases investigated. 

As the solution must be finite everywhere, we therefore have eight 
boundary conditions, namely. six continuity (or jump) conditions at the 
two edges of the plates, and two boundedness conditions at the extreme 
ends of the x-axis. The third-order differential equations applying in the 
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three regions L, M, and N, supply nine integration constants. The damp- 
ing factor, or amplitude ratio, is therefore uniquely determined. 


5. Solution for zero perturbation pressure 

As was mentioned in the introduction, some indication of the behaviour 
of the system can be obtained simply by assuming zero perturbation 
pressure. This would seem to correspond to the case of closely distributed 
waves passing over a relatively long resistance. 

The equation of motion for the uniform basic flow is 


Sian FA, que, (5) 
p 


the bar over the II, indicating that the basic flow is to be taken. Lineariza- 
tion of the first equation of motion for the perturbed flow gives 
n+ , . 
_ Btls _yn(ut42Uu) 
p 


I] 


and therefore Uu, = ——?—hUu. 


If we now assume I], to be of higher order than the first this equation 
simplifies to 
" u, = —hu. 


Therefore the damping factor over the length of the plates is e-*"¢. 


6. Damping factor in the general case 


The general solution of equation (3), valid for |X| < A, is 


3 
= > a,c 


8 


where the A’s are the three solutions of 
84 JA2— PA+41) = 0, (6) 


while the a's are constants. 
For the two outside regions we get, after application of the boundedness 


conditions, , BePXty, (X <—A) 
Bae? +y, (X > A) 
in which the f’s and y’s are constants. 
Substituting these expressions in the continuity conditions at X = +A 


yields a set of six equations, from which the damping factor y,/y, can 
immediately be obtained in determinantal form, as 
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| e-PA —] —e-4 —e-4 ers | 
Pe-PA 0 A, e-A | 
0 —)¢ -\, A | 
_ eA 


r,e4 


AR A\A 


A, A 


—e€ 


A, eAiA 
P2e-PA 0 AZ eA 
in which the blank spaces are filled by cyclic permutation of the index 
of A in the last three columns. After simplification this becomes 
A,(P—A,)e“™"4 A,(P—A,)e“%4 — Ag( P—A,)e-2*4 | 
> 1) 
p?—)i 
P—), 
P?—)? 
—i,(P4 


and since, by (6), A, +A, 


is clearly e~+. Therefore 


4 As(P—A,)(P+A_)( P+Ag)(Ag—Ag)e-™"4 4 
Ax | A,)(P 3 Az)( P -Ag)(Ag —),)e +244 


A form which, for some purposes, is slightly more convenient, is obtained 
by noting that the left-hand side of (6) may be written as 


(A—A,)(A—A_)(A—Ag) = A8+- 4A2— P2(A+1) 


so that 


(P—A,)(P—A,)(P—A3) = —}P? = —(P+A,)(P-+A_)(P-+As) 


and therefore, finally 


V2 = ond A,{(A,— P) (A, 4 P)}(Ag—Aq)e- 2-4 >. 
Y1 Ay{(Ay + P)/(Ay— P)}(Ag—Ag et +...” 
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7. Asymptotic solution for small-scale perturbations 

For perturbation wavelengths small compared with the characteristic 
value 1/h (corresponding to the length of the plates for unit resistance 
coefficient), the frequency parameter P is large, and we obtain the follow- 
ing asymptotic values for the solutions of (6) 


1 
A,= —1 —s + 


Az 


As = 


Substitution into our basic formula (7) gives 


28 ws 4Agl(As— P)/(Ag + P)}A,—Ae ~_ = e4/14+-0( P-? 
é Aol, + P)/(Ay— P)}(Ay —Ag)e244 = e*4/1+4 O(F )}. (8) 





Y1 
8. Solution for large-scale perturbations 
For small P, i.e. for long perturbation wavelengths, the solutions of (6) 
become i, ae hee SOM S..... 
A, = Vv2P— P?+..., 
Ag = —v2P—P?+... 
and substitution in (7) now yields 
v2 Ple4 —4{(V2—1)/(v¥2+1)} —4{(W24 1)/(W2—1)}] an 1—3e-4 


e-A 





V2 Ple-4—4{(v24-1)/(V¥2—1I)}J—H{(W2—1)/(W24-1)}] e438 

(9) 
This expression decreases from 1 to —} as A increases, and vanishes when 
A = \n3 = 1-1. It is valid if A is of lower order than P=}. 


9. Application to distributed resistances and gauzes 
For the interpretation of the results we return to physical variables. 
Integrating (5) over the length of the plates from —a to a gives 
oh we A 
dp? 
which is the drag coefficient k defined by Taylor and Batchelor (3), and 


therefore 2A = Gah = bk 


If / is the wavelength (in the y-direction) of the perturbation, 
p—2"_ 4 


th lk” 
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This means that large values of P correspond to a/k large compared with /, 
representing a widely distributed resistance. The limiting result (8) agrees 
with that obtained by assuming zero perturbation pressure (section 5). 
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If k is kept fixed it is seen that for small P, a is small compared with /, 





which represents a concentrated resistance. The solution (9) for small P 
therefore approximates to the case of a gauze with wires running in one 
direction only. It seems unlikely, when the flow is nearly nermal to the 
gauze, that wires running across the flow would have a constraining effect 
apart from that due to their resistance. The above results for a ‘one- 
dimensional’ gauze may therefore be compared with the results for ‘two- 
dimensional’ gauzes obtained in (3). The resistance coefficient giving 
perfect damping in the present case is k = 2A = 2-2. 

The value obtained by Taylor and Batchelor using empirical relations, 
is k = 2-76. In Fig. 2 the present result is plotted with the experimental 
results of Collar and McPhail (recorded in 3). Over the available range 
of values of k the agreement seems good. It therefore appears that the 
damping effect of gauzes is theoretically explained by the present analysis, 
and that the formulae may with confidence be applied in practice. 
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A NOTE ON LIFTING LINE THEORY 


By K. STEWARTSON (The University, Durham) 
[Received 2 April 1959] 


SUMMARY 
The distribution of circulation round a semi-infinite plane wing of constant chord 
is calculated using the lifting line approximation. The results are used to determine 
the lift and induced drag of a rectangular plane wing of large aspect ratio. 


1. Introduction 

PrRANDTL’s concept of the lifting line has proved of immense value to 
the theory of lift on three-dimensional quasi-planar wings placed at a 
small angle of incidence in a stream of incompressible fluid. In his theory 
the wing is replaced by a straight line /, parallel to its leading edge, around 
which the circulation ['(y) is assumed to be a function of the distance y 
along / from some convenient datum point. Since the circulation is vary- 


ing there will be in addition a trailing vortex sheet, lying behind 7 and 


approximately in the plane of flow of the undisturbed stream, whose 
strength can be determined in terms of ['(y). This trailing vortex sheet 
detiects the main stream in the neighbourhood of the wing, altering its 
effective angle of incidence. On equating this effective angle of incidence 
to the circulation round / by means of the Joukowski hypothesis an integro- 
differential equation for I'(y) is obtained. Once ['(y) is known it is a com- 
paratively simple matter to infer the lift and induced drag of the wing. 

The methods proposed for solving the equation for ['(y) depend either 
on replacing ['(y) by a Fourier series with period equal to the span of the 
wing or on taking a particular form for ['(y) and determining the corre- 
sponding shape of the wing. The reader is referred to Glauert (1) and to 
Robinson and Laurmann (2) for a full discussion of these methods. The 
disadvantage of the methods using Fourier series is that they inevitably 
become less accurate as the aspect ratio of the wing increases while the 
alternative methods have led to solutions relevant to wings with rounded 
tips only. Accordingly it is of interest to determine the form of the 
circulation function for a rectangular wing of very large aspect ratio as 
an example of a wing for which neither of the general methods mentioned 
above is strictly applicable. 

In this note we shall first determine the exact solution for a semi- 
infinite wing of constant chord and then use it to find the leading terms 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960) 
5002 49 E 





50 K. STEWARTSON 


of the expansion of the lift and drag on a rectangular wing, in descending 
powers of the aspect ratio. 


2. The semi-infinite wing of constant chord 
According to Robinson and Laurmann (2), the integro-differential equa- 
tion for the circulation ['(y) around the y-axis is 


’ 8 \ 
: ‘ l I’(y,) dy 
r(y) = Jaye(t Xy + — ( Ty) an, . 

i 4n JO yy 

0 

In this equation the y-axis is in the plane of the wing and parallel to the 
leading edge, the wing extends from y = 0 to y = 8, c(y) is the chord of 
the wing, ay is the geometrical angle of attack of the wing, U is the velocity 
of main stream, and a, is a constant equal to 27 for a flat plate. Further, 

the circulation must be zero at the wing tips and so 


ro) = T'(s) = 0. 


In this section we shall obtain a complete formal solution of (2.1) in the 
particular case of a semi-infinite wing of constant chord so that s = co 
and c is a constant. 

At points distant from the tip the circulation must tend to its two- 
dimensional value and accordingly the appropriate boundary conditions 
on I’ in our particular problem are 


(0) = Q, I(x) 4 daycaU. (2.3) 


Write T(y) = jaycUagfl—f(x)}, where y = aycx/8; (2.4) 


l ‘(x,) dx 
then f(x) = - fan (x > 0), 
7) %%4—2z 
0 
subject to f(0) = 1, f(~) = 0. In order to solve this equation we shall 
first solve formally the slightly more complicated integral equation 


f(a, x) —= | da, f(a, 2,)K,(a\v,—ax!)sgn(z7,—x) (x >), (2.6 


. 
0 


subject to f(a,0) = 1 and f(a,00) = 0, where a is real and positive, and K, 
is the Bessel function of order one, of the second kind and with imaginary 
argument. Having solved (2.6) formally the solution of (2.5) foilows on 
letting a > 0+ and writing 


f(x) = f(0, x). 
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If we set f(a,2) = 0 for all z < 0 it follows that 


F,(w) = [ fla, ae dx 


is a regular function of w in the upper half plane, imw > 0. This is 
expressed notationally by the suffix +. In a similar way a suffix — will 
be used below to signify that a function is regular in the lower half plane, 
imw <0. Inverting the Fourier transform (2.8) 


f(a,x) = us | Powe “tw dey; (2.9) 
2m 


since F,(w) is regular in the upper half plane of w, the integral vanishes 
if~ <0. Substituting (2.9) into (2.6) we get 


a) 


F.(w)+M_(w) = be {—1—iwF,(w)} | dx eK (a x )sgna (2.10) 
7 


— x 

w : 
PN 
(w?+-a?)t u 


= wF,(w)} (2.11) 
where M_(w), like F,(w), is at present unknown, and (w*+-«*)! is regular 
in the w plane cut along the imaginary axis from ia to ix and from — 10 
to —ia. Equation (2.11) may now be rewritten as 

2 


iw 
ae Te, a TT 2.12 
Following the Wiener—Hopf technique, the coefficient of F,(w) in (2.12) 
will now be written as the ratio of two functions, one regular and non- 
zero in the upper half plane imw > —a, the other regular and non-zero 
in the lower half plane imw < «. Thus 
oa ne (2.13) 
L(x, w)  (w* 4-12)? 
These functions may be found by taking a contour (' in the ¢-plane con- 
sisting of the two infinite straight lines C,, C_ parallel to the real axis, 
of which one is just above and one just below the real axis, both lying 
inside jim {| = a. Then if w is a point inside C, 
1 14-72(¢24 
log L.(a,w)—log L_(a,w) = — [ log ia 8 
“7 


— 


— 
. 


Cc 





Hence log L(a,w) = .- [ tog x 
. mi 


- 
7 
’ 
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since this integral still defines a regular function if the range of w is 
extended over the upper half plane. Differentiating (2.15) and integrating 
by parts we find 

(6/6w) L(x, w) . l i 


2.16 
L (a, w) ( ) 


whence L.(«,w) may be obtained formally and shown to satisfy the 
conditions required of it. However, in the present instance we are most 
interested in the particular case 1 = 0 when the integral reduces to 
1 / Cdl 
m J (C?—w*)(f+1) 
0 


(2.17) 


The two integrals by which L,(0, #) may be found are most easily evaluated 
when w = it, t real and positive, and we then find that 


t 
1 f log é dé| 
a) 1462 | 


(2.18) 


L.(0,it) = (1 HA exp| 


taking L .(0,0) L .(x, 0) 1. Having obtained L ,(0,w) on the positive 

imaginary axis it can easily be evaluated elsewhere. It is only necessary 

to cut the w plane along the negative imaginary axis. If, for example, w 

1S real ; w 

L (0, w) (l+- w ree | log a (2.19) 
1— 6? | 


77 
. 


0 
Similarly, it may be shown that if w iu, u real and positive, 


“ 


] [ log @ dé| 


2.20) 
142 | \ 


L_(0,w) (1+-u?) texp! 
0 
Having formally determined L.(x,w), L_(a,w) let us now reconsider 
(2.12), which may be rewritten as 


lw 


F.(w)L .(x,w) L_(a,w)—M_(w)L_(a.w). 


(w? t x”)! 
and hence, using (2.13) to eliminate (w?+-«?)~?, as 


F.(w)L.(a,w) : {Lo (a,w)—1]} L_(«,w)M_(w) ‘ {LL _(a,w)—1}. 


w Ww 
(2.21) 
If now we assume that F,(w) is regular in im w > —a and M_(w) is regular 
in imw +a it follows that both sides must be equal to a function H 
of w which is regular everywhere. The form of H can be determined by 
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considering the behaviour of F, and L, as w+. Since f(a,7)—> 1 as 
xr->0, wF, >i as w—>o. Again from the symmetry in the definitions 
of L.(a,w) and L_(a,w) in (2.13), L.(«,w)/wt must tend to a finite limit 
as w-» oo. Hence as w +o the right hand side of (2.21) must tend to 
zero, so that H(w) = 0. The verification of the assumptions stated at the 
beginning of this paragraph now follow a posteriori and we have 


i | t l 
Pr, t—--—) pe Sfp, _) ess 

(w) +( men) () w'\ L all ( ) 
The formal solution of (2.6) is now completed by using (2.9). Of particular 
interest is the special case a = 0 for which L, and L_ have been found 
as (2.19), (2.20). In this case 


x 


] i l iwsr 
fie) = fi0,2) = 5 [Al ia ” 


L t 

lf e# ie 
-exp| — 

; j 


log @ d0) 


» Oo 
1 j #2 | dt, 2.23) 


ape *P| | 
0 
on deforming the contour into the two sides of the negative imaginary 
axis of w and using (2.13), (2.20) to find the values of L,.(0,w) on either 
side of the cut along this line. From (2.23) it is a comparatively simple 
matter to compute the values of f(x) displayed in the table below: 





x oO “2 o'4 06 o8 ro : ° 6 “s . 30 40 





f(x) 0° 564 | 0°438 | 0° 361 | 0°308 | 0267 4 "175 O11 3 | o-086 
4 4 j 












































If x is large, f(2) can be found by expanding the integrand of (2.23) in 
ascending powers of t and log? and integrating term by term. It is found 
that 


[ 2(log x)? +-4y log r—1—$n?]+... (2.24) 


9 
< 


x 


f(x) 


, (y +log x) + 


— 
77 


mx) 
as x » x, where y is Euler’s constant 0-5772.... On the other hand, if x 


is smal), f(x) is found by expanding L ,(0, w)in descending powers of w and 
interpreting each term separately although some care is needed. We 


f(x) =1 2(7)! il.) Mosse ty yea) (2.25) 
\7r 3\r 7 


where 0, means that the order includes an unspecified power of log x. 


obtain 
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Perhaps the most useful property of f is the value of 
[ f(x,) dx, 
when 2 is large. From (2.9), (2.2: 


l 
L (0, it) 


[ste wt dt 


if ¢ is real and positive. Now 


a(x-+-1) 
0 


e “dx Af tlogt)4+ O(tlogt) (ft small), 
7 
and therefore 


y y+] 
| F(a) a0 . ax lim | 3 ‘ 
hs a nm(a-+1) t 


— | 1% ot : g 
1L.(0, it) ta 41 89. 


t--0 7 
0 


Accordingly, when is large, 


1 f( r,) dx, (log x +y-+-1)+-0(1). 
F 7 
Finally, using numerical methods, it was found that 


0-32. 


3. The lift and drag on a rectangular wing of large aspect ratio 


Here we consider a wing of constant chord ¢ and span s such that 
A sc » 1. For this wing write 


( , 
0 


(y) = 
where 8s a,cr. Then 


* g'(x,) dx cE Fa 
g(x) Bt “Anand | | / (x,) da, 
ry P 7 . vr 
r 


and g(0) = g(r) = 0. However, if r is large, f(x) 
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hence the last integral is O,r-*). It follows that g(x) = O,(r-?). The lift 
on the wing is 


L pU | Py) dy = pU*%esC, 


ay C7 xy pl 
* Té6i rat ee fir) f(x) f(r —z)+ g(x)| dx 


A C2x9 pU® f14 — Aloe r+-1-+-v)-+-Ofr- 
- 16} l t-f(r)} slit }r 2( ogr y) ir ). 


1 8A 
Hence C, =a I 40 _llog"“ 4.1 | tO, A-? 3.5 
L 0 0 4nA 2 a y¥| +A ) (3.5) 
when A is large. Taking a, = 27, (3.5) gives a value of 4-99 for the slope 
of the lift coefficient at A = 10 which may be compared with a value 5-04 
quoted by Glauert (1). Equation (3.5) should also be compared with 


Cy = 4% vf 9 4Q(A »| (3.6) 
aay, 


which is the equation for lift coefficient of an elliptically loaded wing of 
large aspect ratio. This means that the effect on the lift of having 
rectangular instead of rounded tips to a wing of large aspect ratio is 
manifested by the presence of the logarithmic terms in (3.5). 
The induced drag of a three-dimensional wing is 
‘ 


[ pul (y) dy = tpU%esCp, (3.7) 


dD, 


(1, 140), where w is the downwash and is connected with the circulation 


through the equation WT = aie—wi. (3.8) 
Substituting (3.1), (3.8) into (3.7) and using (2.30), (2.31) it is found that 


2,2 
C Ay & 


8A ‘ 
D, oe (log = +y4 0-0) + 0,(A-*) (3.9) 


when A is large. Expressing this result in the manner conventional to 


aerofoil theory, 


cy 8A 
Cy, = eet. 0, 1 3.10 
n= Al (loss +7) +004-] (3.10) 
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when A is large, which may be compared with the formula 


for an ellipticaliy loaded wing. 
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A SOLUTION OF THE COMPRESSIBLE LAMINAR 
BOUNDARY LAYER EQUATIONS WITH HEAT 
TRANSFER AND ADVERSE PRESSURE GRADIENT 


By G. POOTS 


(Department of Theoretical Mechanics, Bristol University) 
{Received 12 August 1958; revise received 19 December 1958] 


SUMMARY 


A solution of the two-dimensional compressible laminar boundary layer equa- 
tions with heat transfer and adverse pressure gradient is given. The problem arises 
from the application of Stewartson’s transformation to the compressible laminar 
boundary layer equations. In the transformed incompressible-plane the main stream 
velocity and dimensionless temperature at the surface have been chosen to be 
U, = u,o(1— 4X) and S, 1 respectively, where X is the longitudinal coordinate 
in the incompressible plane. These forms have been chosen so as to simplify the 
numerical work involved in solving the coupled partial differential equations which 
oceur. A solution for the coupled partial differential equations is obtained, for small 
X, by means of a power series in terms of X, where the coefficients in the expansion 
are obtained on integrating a series of ordinary differential equations; for values 
of X approaching the separation point the partial differential equations are inte- 
grated using standard relaxation methods in conjunction with the Hartree— 
Womersley method. 

In the incompressible-plane boundary layer characteristics such as displacement 
thickness, ete., are evaluated from the numerical solution and compared with those 
obtained by considering an approximate solution of the transformed momentum, 
energy, and thermal-integral equations. On transformation of these results to the 
compressible-plane information is obtained for the compressible laminar boundary 
layer on a flat plate with heat addition and in the presence of a retarded main stream 
velocity distribution. In particular dimensionless quantities related to the shear 
stress and heat transfer at the wall are evaluated at various stations along the plate 
from the leading edge to the point of separation. 


1. Introduction 


Tue solution of the compressible laminar boundary layer equations de- 


pends on the pressure gradient in the main stream, Mach number, and 
heat transfer over the solid boundary, together with the properties of the 
fluid under consideration. As shown by Cope and Hartree (1) the numerical 
work involved in a complete study of these equations is prohibitive when 


allowance is made for the empirical variation of viscosity with tempera- 
ture. However, Howarth (2) and Stewartson (3) gave a transformation 


such that the boundary layer equations of the compressible fluid may be 
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transformed into a form identical with those of the incompressible case. 
In Howarth’s analysis the following assumptions are made: 


(a) the surface of the body is insulated, 
(b) the viscosity («) varies as the absolute temperature 7’, 
(c) the Prandtl number (c) of the fluid is unity. 


Stewartson (3) further simplified Howarth’s analysis and extended it 
to include the case of arbitrary surface temperature. It is the purpose 
of the present paper to give a numerical solution of the compressible 
laminar boundary layer equations for the flow along a flat plate, main- 
tained at a certain surface temperature, and in the presence of an adverse 
pressure gradient. In the following we adhere to Stewartson’s analysis. 

The equations of the steady, two-dimensional compressible laminar 
boundary layer for perfect fluids are, in the usual notation, 


ou ou Cp of eu 
pu at 3 pv- » = oa a p—}> (1.1) 
CX cy ce OYy\ cy 


Gu ——, (1. 

cy 
~ (pu)+—(pr) = 0, (1. 
Ox Cy 


and Je, pu: ae pv: *) "oe (Yeon) al: s) (1. 


Cx cy cx Cy oy oy 
These equations may now be transformed into equations similar to those 
for the two-dimensional incompressible laminar boundary layer by means of 


Stewartson’s transformation. Stewartson introduced the following variables 


dX =d x| 


‘A 


dY =—dy * (2); 


Ag \V¥o\Po, 


and the stream function defined by 


=” I)(y-1) 


ct, §, wh, (1.7) 
Po Yo Po *¥o 
where the suffix 0 refers to some standard state, taken to be the main 
stream at x = 0. The suffix e refers to local conditions at the outer edge 
of the boundary layer. Applying equations (1.5) to (1.7) to the boundary 
layer equations (1.1) to (1.4) and assuming that o = 1, we obtain the 
following equations in the incompressible-plane 
eu. av 


CoS — 0, 1.8 
ox * oY ati 
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éU  .,aU dU, 
Ux toy = U-=¥ (I +8)4+- 5%. 


yes , yes _ #s 
aX" @Y ay?’ 
where S is the dimensionless temperature defined by 


T Ti (y—)) yalg. Y¥—)) ‘ 

T= +142 MoS + saa (ue—u*), (1.11) 
and where the stream function has been replaced by the transformed 
velocities (U, V) defined by 


U=y%, V=—d,. (1.12) 


The resulting relationship between the transformed and physical velocities 
is U = (a,/a,)u. In the main stream we have also that 


“4 of 4 N) 2 = a3 + %— Mya = constant (1.13) 


— 


and (" —_ = (1.14) 
Po 4 

The relationship between U, and u, is [1,1], so that for every problem 
in the incompressible-plane there is an associated problem in the com- 
pressible-plane. Thus we have two alternatives: either we can take an 
assumed form for the retarded main stream velocity u, in the compressible- 
plane and determine U, in the incompressible-plane by equations (1.5) and 
(1.13) or vice versa. Since the assumed form could never be obtained 
physically, it is advantageous from the viewpoint of the numerical work 
to assume a simple form for the retarded main stream velocity in the 
incompressible-plane. Hence we shall consider a transformed main stream 
velocity distribution given by 


U, = uo(1—}X). (1.15) 


Furthermore we shall assume the dimensionless temperature at the wall 
in the incompressible-plane to be given by S,, = 1. The main stream 
velocity and temperature distribution, together with the wall temperature 
distribution in the compressible-plane may now be obtained. Substituting 
equation (1.15) in (1.13) we find that 

. ij4 —1)M? 
Ge\* _. tH DM (1.16) 
at) 1+hy—-I)M, 
where M, = M,(1—}X). Therefore 


le ey , 
we = “¢(1—4X) (1.17) 


0 % 
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4 2) (8y—1)/Xy-1) 
and z= "ltho—1) mol pies dX 


yD 


ake 


If y = 1-4, the integral (1.18) reduces to 
aw = (14+}M})-4[X+2MF1—(1—}X)3}+ BMe{1—(1—3.X)}+ 
+ 3% MS{1—(1—3.X)"}+ ghsg ME{1—(1—}$.X)*}]. (1.19) 


In the main stream S 0; hence 


7S 14 (Y- D (ud —ud 
T = 2a2 e 0 
e ~~ 


), 


giving, on using equation (1.16), 
T T,- a i a (1.20) 
+4(y—1) MZ 
The dimensionless temperature distribution S evaluated at the wall deter- 
mines the wall temperature. Thus 
T= 1+ gh HD MBS + Hy ME, 


and on using equation (1.20) the result is 
T,, T(1+S,){1+4(y— 1) Mj}. (1.21) 


Hence the case for which M, = 1, y = 1-4, S,, = 1 corresponds to a wall 
temperature 2-4 times the free stream temperature 7), at x = 0. 
The boundary conditions to be imposed on equations (1.8), (1.9), and 


(1.10) are U(X, 0) = V(X,0) = 90, 
S,, = S(X,0) =1, 


w 


lim U = U(X) = u,(i—4X), 
Yy—« 
and lim S 0. 
Yau 


Thus by solving equations (1.8), (1.9), (1.10) subject to (1.22) information 
on shear stress, heat transfer, etc., may be obtained for the compressible 
boundary layer on a flat plate maintained at wall temperature 


T,, = 2T{1+}My—1)MG} 
and in the presence of a retarded main stream velocity u, = uo(1—}X)a, ay 
and main stream temperature distribution T, = 7),(a,/ay)?._ In particular, 
separation will occur in the compressible-plane when (Cu/éy),<9 = 9, i.e. 
in the incompressible-plane, when (éU’/é¥)y.5 = 0, since the shear stress 
evaluated at the wall 


eu a,\4y-2ky-D (aU ia 
Tw (1 Po YY = . (1.23) 
Cy] y-0 Ay OY} y <0 
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In section 2 of this paper a series solution in powers of X for the coupled 
partial differential equations (1.8) to (1.10), subject to the boundary 
conditions (1.22), is obtained following the method used by Howarth (4) 
in his discussion of the problem of the incompressible flow over a flat plate 
in the presence of a linear retarded main stream velocity. The coefficients in 
the expansion are obtained on the numerical solution of a linear set of 
ordinary differential equations. The amount of numerical work involved 
is large and the expansions have only been evaluated up to and including 
the term X*. This is insufficient to obtain precise information near the 
separation point. In section 3 details are given on the determination of 
velocity and thermal profiles at various values of X, up to the separation 
point, by solving the basic equations using a numerical method developed 
by Hartree and Womersley (see Hartree (5)). Section 4 deals with the 
approximate solution of the momentum, energy, and thermal integral 
equations in the incompressible-plane by means of a modified Polhausen 
method. The boundary layer characteristics, such as displacement thick- 
ness, etc., so obtained are then compared with those evaluated using the 
more precise results obtained in sections 2 and 3. In section 5 the dimen- 
sionless quantities u/u,, T,/T), T,,/T, and those related to the local skin 
friction C, = 7,,/$p,,u? and the wall Nusselt number 


oT 
N — T.—T.)-1 
y (7), s “ w) 


are calculated for the compressible-plane for a range of free stream Mach 
numbers (J, = 1, 2, 3, 4, and 6) as functions of x/x,, where x, is the 
distance of the point of separation from the leading edge. Finally in 
section 6 conclusions are given. 


2. Series solution 

We consider here the solution in series for the boundary layer equations 
in the incompressible-plane. Following Howarth (4) we introduce the 
variables é= X, 7 huh YX. 2.1) 


The stream function is defined as 
b — up Ed(€, n), 
where the function ¢(€, 7) must satisfy the equation 
git, of BAY ob 
3 2T 
en oe c& 7? ~ én e&On 
and the dimensionless eines S must satisfy 


28 ‘ ‘ 
ae -24(4 28 298) _ a 
cE On On CE 


)— —£(1—4£)(1-48) = 0, 
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The boundary conditions are 


d op 0, Ss 1 at » 0 for € >0 


C7) 


lim % —2-3¢, limS =0 foré > 


t/a c 7) aks 
We now assume expansions of the type 
$(f,4) = > (—1)'€9,(n) 


r=0 
and S(E, n) > (—1y ES»). 


r=0 


When equations (2.6) and (2.7) are substituted into equations (2.3) and 
(2.4) and like powers of € are equated, the following differential equations 
result: 

$9 +b5% = %, (2.8) 


$” +-45¢1—2654,+3¢04, (1-4.8,), (2.9) 
$5 + bob5—445 43+ 585 by — —M1+S,)—S,+2(¢1)?—34, 4%, (2.10) 
and for r 


d, + do d, -2rdo d, i-(2r- 1)d6 4, 


x B 
ra Sat SX (246, ¢)—(2P+ Vty $i}, (2-11) 


pe-lgq 
where the prime denotes differentiation with respect to 7 and the summa- 
tion is restricted by the selection rule a+8 = r and a, 8 = 1, 2, 3 
Equation (2.4) for S reduces to 


r. 


gees 


x Fe] 
2 t (2r H l dy S, 2rdyS, > 2! 294,58 (2p 1)¢, S3 (2.12) 
p=1q 
for r = 0, 1, 2,..., where the double summation is restricted by the same 
selection rule as its counterpart in equation (2.11). The boundary condi- 
tions are 
for r 0, 1, 2,... at » = 0, 
for r 1, 2, 3,... at » = 0, 
i, d, 0 forr 2, 3, 4,... at 7 00, 
and i for r 0, 1, 2,... at 7 00. (2.13) 


The differential equations for ¢,, $2, 43, ete., and Sy, S,, S,, ete., are 
ordinary linear non-homogeneous differential equations with variable 
coefficients. Since the coefficients on the left-hand side of each equation 
are not known analytically, the ee must be integrated numerically. 
The functions ¢,, 45, 43, $4, S,, S,, S,, and S, and their derivatives have 
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been obtained using the Adams—Bashforth method, and the Blasius func- 
tion ¢, and its derivatives are as tabulated by Howarth (4). The calcula- 
tions were performed to seven decimals and the resulting functions quoted 
to four decimals are given in Tables 1 and 2. The function S, is not 
tabulated as it may readily be obtained from Howarth’s table of 49, since 
So 4(2 do). 


In the incompressible plane separation will occur at X — X, when 


v “9 
2 
7 4) no 
on 7=9 
Although all the ¢$7(0) increase with r in absolute magnitude, some 
indication of the value of X, may be obtained from the existing numerical 
work. Separation occurs when 


a 

> (—1X"45(0) = 0. (2.14) 

r=0 
Noting that the ratios |43(0)! //43(0)| and /43(0)| / 45(0)! are approxi- 
mately equal to 1-02, let us suppose that |47(0)| = 1-02'47_,(0)| for r > 4. 
On using Table 1 a suitable approximation to equation (2.14) is then 


4 
1-32824— 1-5939X — 0-2608X —0-2644X2—0-2780 =0, (2.15) 


giving an estimate of X, = 0-64. The above approximation to equation 
(2.14) presupposes from the slender evidence already obtained that the 
',(0)| do in fact increase as r increases. It is interesting to note that 
more detailed numerical work, as carried out in section 3, gave X, = 0-60. 
This indicates that the '47(0)| may increase more rapidly for large r than 
estimated above. 

Although the above series expansions (2.6) and (2.7) could, with some 
labour, be extended beyond the term X*, it was found that considerable 
accuracy was lost in proceeding to the calculation of the higher functions. 
This was mainly due to the numerical magnitude of the variable coefficients 
on the left-hand side and to canceilation within the summation on the 
right-hand side of equations (2.11) and (2.12). Nevertheless the finite 
forms 


=> (-1V€'4, 


r=0 


and S = > (—1)€8, (2.17) 
r=0 

were found extremely useful in checking at € = 0-2 the velocity and 

thermal profiles obtained by the Hartree-Womersley method, as described 

in section 3. 
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3. The application of the Hartree-Womersley method to the basic 

equations (2.3) and (2.4) 

The essential idea of the application of the Hartree-Womersley method 
(see Hartree (5)) to the basic equations (2.3) and (2.4) is that at any two 
stations € = &€, and € = &, we replace derivatives in the é-direction by 
finite differences (without truncation error) and all other quantities by 
averages. Thus we obtain 


($4 +p) + ba{(b.4)?—(bp)*}+ O44 belbat+bnt bp—by)}+ 
+a(24+5,4+S_\URe—U%e) = 0 (3.1) 
and 
(5% T Sp) T ha(Sy Subs +p) t ba( 8’, + Sr)(op 3 4) T 
t¥(S'4+Sp\(di+bn) = 9, (3.2) 


where « = 2(€,+&,) (€,—€&,) and the primes denote »-derivatives at the 


stations € — €, and € = €,. The boundary conditions for ¢,, dz, Sy, 
and S, are 


oy dp = ¢, d'; t ’ i atyn=—0 foré>0, 


2— hE 4, dy = 2—}ép, S S,=0 atyn 0 for € > 0. 
(3.3) 


If we now introduce 


S=¢+¢, ¥ , 0=8,48, (3.4) 


then equations (3.1), (3.2), and (3.3) simplify as follows: 
W"+-{4(1+a)D—ad,}" +a(b)y —§V)¥ +-a( U§,e—U*, e)(24-6) = 0, (3.5) 


, 


8” +-{3(1-+-a)D—ad,}0’—Ja¥0+a5,¥ = 0, (3.6) 


together with the boundary conditions 


7 = 6, 6=2, aty= 0, 

Y = 2(U,e—U,e) = 4—}(Eqt+Es), 9= 9 at y = @. (3.7) 
Equations (3.5) and (3.6) represent a pair of coupled ordinary differential 
equations in VW and @, the equation in ‘VY being non-linear. If we know ¢,, 
¢4, and S, the boundary value problem is then completely defined. 

The integration procedure is started at €, = 0 since ¢,, $4, and S, are 
known from Howarth’s table (4). Taking €, = 0-2, and thus a = 2, the 
required functions ‘W and @ may be found by standard relaxation methods. 
These functions VY and @ will only be approximate due to the use of simple 
truncated finite difference formulae in the €-direction. An improvement 
in accuracy is then obtained using Richardson’s process of h?-extrapola- 
tion. We calculate ‘VY and @ at €&, = 0-2 in one and two steps € = 0 to 
0-l and € = 0-1 and 0-2. If WY and @ are the results for the one-step 
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calculation, and if ‘¥Y@ and @ are the results of the two-step calculation, 
then by A®-extrapolation the ‘final’ results are 


Yy — PO Yea yo) 
and 0 = G)4-4(G2)— A) (3.8) 


with error O(h*), where A is the mesh length in the €-direction (see Fox (6)). 
That h?-extrapolation is the correct extrapolation for the above equations 
has been verified numerically by comparing at € = 0-2 the functions ¢, ¢’, 
and S (obtained from the final h?-extrapolated functions ®, ‘Y’, and @) 
with the series solution given by equations (2.16) and (2.17). The difference 
was found to be not more than one unit in the fourth decimal, which is 
in good agreement with the allowed possible error O(/*), since h == 0-2. 
Further details on the relaxation solution of equations (3.5) and (3.6) for 
& < 0-2 and 0-2 < € < &,., is described in an appendix. 

In Table 3 the velocity U = }u,é¢/én, and thermal profiles S are 
tabulated at intervals of 7 = 0-1 for the stations € = 0-2, 0-4, and 0-6, 
At stations € = 0-2 and 0-4 the functions U’/u,g and S are rounded off to 
four decimals since the difference between the one- and two-step calcula- 
tions was less than ten units in the fourth decimal, the maximum difference 
occurring at 7 = 1-4. It is believed that the fourth decimal in the ‘final’ 
h?-extrapolated values of U/uy and S, at € = 0-2 and 0-4, is reliable. 
However, in the calculation of U/u, and S at € = 0-6, in one and two steps 
from & = 0-4, the difference in the one- and two-step calculation was less 
than forty units in the fourth decimal, the maximum again occurring at 
about » = 1-4. Final h?-extrapolated functions at € = 0-6 are thus rounded 
off to three decimals and it is believed that the third decimal is reliable. 

At & — 0-6 the shear stress, (@U/cY), 9, in the incompressible plane was 
found to be of the order 10-4. Hence separation will oceur very near to 
X, = 0-60. Since the calculations at € = 0-6 become rather unstable, it 
is possible that the numerical process may break down if an attempt is 
made to pass the point of separation. Indeed, for the incompressible flow 
over a flat plate in the presence of a linear retarded main stream velocity, 
Hartree (5) found, on reinvestigation of Howarth’s problem (4), a strong 
suggestion of the existence of an algebraic singularity in the neighbour- 
hood of separation. This conclusion was later confirmed by Leigh (7) in 
repeating Hartree’s calculations using much finer steps in the €-direction 


(in the neighbourhood of the separation point) and working to a greater 
number of decimals using automatic computing machinery. The question 


as to whether or not an algebraic singularity in (¢U//éY)y_» exists for the 

present problem remains unanswered due to the prohibitive amount of 

necessary numerical work. It has been the intention in this section to 
5002.49 F 
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obtain over-all trends and gross effects rather than precise results. How- 
ever, the velocity and thermal profiles obtained are sufficiently accurate 
to examine the exactness of various approximate methods which may be 
developed and which proved so fruitful in the study of the incompressible 
boundary layer. 


4. The momentum, energy, and thermal-integral equations in the 
incompressible-plane 

In the previous two sections details have been given on the numerical 
solution of the transformed boundary layer equations (1.8), (1.9), (1.10), 
and (1.20). It is desirable to compare with these solutions the accuracy 
of approximate methods. Since any approximate method applied in the 
incompressible-plane can be taken over directly to the compressible-plane 
(see Stewartson (3)) we shall develop, for the sake of simplicity, an approxi- 
mate method in the incompressible-plane. For the laminar incompressible 
boundary layer the skin friction and dissipation of energy are connected 
with the boundary layer thickness by the momentum and energy-integral 
equations, which represent the balance of momentum and energy within 
a small section of the boundary layer. These equations are due to von 
Karman (8) and Wieghardt (9). In the following we shall derive the 
momentum, energy, and thermal-integral equations in the incompressible- 
plane and proceed to obtain an approximate solution of the transformed 
boundary layer equations (1.8), (1.9), (1.10), and (1.20). 

We make the assumption that, in the incompressible-plane, the hydro- 
dynamic and thermal boundary layer thicknesses are identical and equal 
to A. Integrating equations (1.9) and (1.10) with respect to Y between 
the limits Y O and Y A, and taking account of equations (1.8) and 
(1.20), we obtain the momentum and thermal-integral equations 

7 49" , 2. galo . (A*- E)\ F Tr ) 
“dX dX | Oo | U,\e@Y] y-6 


e 


(4.1) 
a” aa eY 


Multiplying equation (1.9) by U and integrating with respect to Y between 
the limits Y = 0 and A, we have the energy-integral equation 


A 
,do*? |, dU, 2E*)  40* f (aU? 
l . » he 2B € (=) #2 * —_ f 1y : 43 

* Cn dX | ae) v2 } (“) C ( ) 
( 


” dE*  @U, Re 2E*() (4.: 
¥y-0 


Here the boundary layer characteristics in the incompressible-plane, i.e. 
the displacement thickness A*, momentum thickness ©, energy thickness 
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©*, thermal thickness Z, and thermal flux thickness £*, are defined as 


A 
At — ( -7) ay, (4.4) 


€ 


Vid su ™ 
ae (4.5) 


€ 


(4.6) 


4 
.{ say, (4.7) 
5 


E* SdY. (4.8) 


r 
. ( € 
( 
Note that if S is put equal to zero in equations (4.1) and (4.3) we obtain 
the momentum and energy integral equations derived by von Karman (8) 
and Wieghardt (9) respectively for the incompressible boundary layer. 
In the original method of von Karman and Polhausen (10), applied to 
the incompressible boundary layer, a fourth degree polynomial for the 
velocity profile was assumed and the coefficients were determined using 
the boundary conditions at the wall and in the main stream. The momen- 
tum-integral equation could then be solved approximately. It was found 
that the velocity profile could be expressed as a member of a one-para- 
meter family of curves, with A = ©?(dU,/dX), the so-called Polhausen 
parameter. On examination of the various exact solutions known for 
incompressible boundary layer theory, it was found that the above pro- 
cedure was adequate for favourable pressure gradients but for adverse 
pressure gradients it is evident that the Polhausen parameter is not suffi- 
cient to determine the velocity profile. Wieghardt (9) suggested using a 
two-parameter velocity profile, one parameter being the Polhausen para- 
meter, to be found in the usual manner from the momentum-integral 
equations, and a second new parameter which was to be found by using 
the energy-integral equation. In Wieghardt’s method an eleventh-degree 
polynomial for the velocity profile was assumed. This procedure, although 
it gives satisfactory results, has been further simplified by Tani (11) with 
considerable improvement in accuracy. Tani assumes a fourth-degree 
polynomial in which four of the coefficients are determined using the 
conditions at the wall and in the main stream and the remaining coefficient 
is determined by satisfying the energy-integral equation. In what follows, 
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Tani’s method is modified to deal with problems in the incompressible- 
plane. 
We assume velocity and thermal profiles of the type 
4 4 
Bl YX 7 Y 4 ~ = ¥ n 
> a,(X)i—] . S - b(X 
=D, @atX)(5] > Pal 
n=0 n=0 
and determine the coefficients a, and b,, from the conditions 
U=0, S=-1 a&@& Y2xO, 
U C 21 ] cS es 
U=U, Ss=0, : -=0 atY=A, (4.10 
; oY oY? oeY eyY? \ 
The usual conditions 
2U dl’ ,dl, (78 
: ") ut €(1-4-S8(0)) out ¢, athe _ 
oY? Y=0 dX dX ey? ¥Y=0 
which imply that equations (1.9) and (1.10) are satisfied at the wall, are 
not used so that in each case one of the coefficients denoted by a and 6 
respectively may be left undetermined. We thus obtain 
U "2 P 2 a "\3 
‘(5-8 3a) +e5 1-5] (4.11) 
U, AN A A’, A A 


and S ] -7(6 an y 1.3 us _»? _ I ). 
A? A A? A A 


The ‘new’ coefficients a and 6 are then adopted as parameters for the 
velocity and thermal profiles, and are to be identified with the velocity 
and thermal gradients at the wall since 


re 


ev S| 
: a =. (: (4.13) 
U\e¥]y.9 =A eY]y~0 
Introducing the dimensionless quantities 
-) @* 
€= ' e f= 
A A 
A 
2A [ (' Ae ao 
‘Y, 
U? | C 7) 


? 


E* 
and B P 
(-) 


we can write egnations (4.1), (4.2), and (4.3) in the forms 


l ail { 20°F (2 i q., “| = p, 
“dX dX \ e ¢€ 


ye, 2,202 (3 +2") = 9, 
dX dX \ f 

dU, _ 

dX 


I d(B*)*) | 
dX — 


and - 28°02 r. 
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The quantities d, e, f, g, h, k, p, q, and r are functions of a and 6 only. 
Substituting equations (4.11) and (4.12) in equations (4.4) to (4.8) and 
using the definitions given in equation (4.14), there resultst 


d 2 @ 
o= -——, 
5 20 


4 a 


35105 


876 


ae A Ee 
5005 | 5005 5460 


2 »b 
cages 
5 20 
4 25 


35. 840 


a® 


cates 
252 


5D l 


=a 2 , _ qi 
2860 


17 


840 


9 
: .~ (48—4a + 3a?), 
35 


€ 


2 (aU . 
— = 28 ’ 
at r), Ba 


(4.18) 


(4.19) 


(4.20) 


(4.21) 


(4.22) 


(4.23) 


(4.24) 


A 


{72 


€ « 


4-)* ( U 


, ) ay kf, (4.25) 
a 


(eS ; 
and = —2 r( ) == 2bh. (4.26) 
EY] yo 
The first-order differential equations (4.15), (4.16), and (4.17) are to 
be solved subject to the boundary conditions © = 0 at X — 0. They 
represent, together with equations (4.18) to (4.26), three coupled first- 
order differential equations in a, 6, and A. A more suitable form of these 
equations for numerical integration is obtained on introducing the Pol- 
hausen type parameter 
dU’ 
n2A = ai? —<£. 
dX 
Then for the main stream velocity U, = u)(l1— X), the approximate 
energy-integral equation (4.16) becomes 


. iA h 
x .@—+Al6+451 = @ 2) 
a % (6 7) q (4.58) 


t+ Formulae (4.18) to (4.20), (4.23) to (4.25) have already been given by Tani (11). 
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and on using equation (4.16) we can eliminate d@?/dX from the equations 
(4.15) and (4.16), giving 


ke ef h (d+g)\ =. ) 
ee ae WN 4s A(s—X)— — [log 
; f fi f é | (s 3 "ee 


and b= kh (1 h 


f° f? pee 


ff} 
Since a = ©* @ is finite at X — 0, the boundary condition to be imposed 
on equations (4.28), (4.29), and (4.30) is 


d - f , 
ax(!°8 5): os 


A=0 at X¥=60, (4.31) 
The integration of equations (4.28), (4.29), and (4.30) subject to equation 
(4.31) has been carried out using the Milne-Simpson step-by-step integra- 
tion procedure for first-order differential equations (see Milne (12)). The 
interval was taken to be X = 0-02 and the iterative procedure for calcu- 
lating a, }, and A was carried out using four significant figures. Although 
it is difficult to extend the calculation right up to the separation point 
(i.e. when (cl CY),-~ aA 0), due to the rapid growth of the deriva- 
tives d(logf «) dX and d(logf h) dX in this region, it is possible to arrive 
at a point so near separation that the latter point may be extrapolated 
with confidence. It was found on extrapolation that 


feu’ a 


| | 0-000 at X 0-60, 
6Y}/y-9 A 


i.e. a O(10-%), 


Using the tabulated values of a, 6b, and A at X 0(0-02)0-58 and the 
extrapolated values at X -- 0-60 we can calculate, for the incompressible- 
plane, the boundary layer characteristics 


A* 6, 6*, Zz, £*, (5) and (sy . 
ey Y-0 Ug ey Yy=0 


In Table 4 these approximate results are tabulated for the stations 
X = 0-2, 0-4, and 0-6 together with the more precise values evaluated 
from Tables 1, 2, and 3. We note that the agreement is excellent at 
X = 0-2 and 0-4 and the discrepancies at X = 0-6 are not serious. Results 
for the approximate solution are given graphically in Figs. 1 and 2, and 
those obtained using the series solution and the Hartree~-Womersley 
method are the indicated points in these figures. 


5. Boundary layer characteristics in the compressible-plane 

In sections 2, 3, and 4 the boundary-layer characteristics in the in- 
compressible-plane have been calculated for the main stream velocity 
U, uo(1—}X) and dimensionless wall temperature S,, = 1. It is the 
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Fic. 1. Hydrodynamic boundary layer characteristics in the 
incompressible-plane. 
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2. Thermal boundary layer characteristics in the incompressible-plane. 
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purpose of this section to give in tabular and graphical form some of the 
essential boundary layer characteristics in the compressible-plane for a 
range of free stream Mach numbers M, = 0, 1, 2, 3, 4, and 6. The argument 
of tabulation is taken for convenience to be x/x,, where x, is the distance 
of the separation point from the leading edge. In the following we take 
y ]-4. 
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Fic, 3. Main stream velocity distributions and the associated main stream 
temperature distributions for My = 9, 1, 2, 3, 4, and 6. 


(i) Velocity and thermal profiles in the main stream 

These have already been defined by equations (1.17) and (1.20). The 
quantities u,/u, and 7./T, are given graphically in Fig. 3 and in tabular 
form in Table 6. 


(ii) Shear stress and local skin friction coefficient 
The shear stress at the wall is defined by equation (1.23). On using 
(2.1) and (2.2) we obtain 


% (52, 
= us (7d - 
Tw = Po \ooaxilena) . (5.1) 
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This can be made dimensionless on introducing the local skin friction 
coefficient C, = 7,,/4p,,u2 and the Reynolds number Re,, = xu,/v,, result- 
ing in the relationship 


' , \B8y-—DAAy—D} /A2, 
Cyy(Beeg) = dettx— pays a(S ES) 
7) én*) »-0 


1 +(y-1)M,? Nu 
1+%*M,? VRew 








x/; Xs 
> 








0001 02 os 04 05 068 OF 08 09 
Fic. 4. Effect of an adverse pressure gradient on wall shear and 
Nusselt number for 0 < My < 6. 


The negligible variation of C,,(Re,,) with Mach number will be noticed 
from Table 6. For this reason the single curve of C,,(Re,,) against x/x, 
in Fig. 4 is sufficient for the range of Mach numbers 0 < My, < 6. 

(iii) Nusselt number 


The local heat transfer at the wall is described by means of the dimen- 
. , * oT ry a 
sionless Nusselt number Nu o(' jin T,,.). The Nusselt number can 
ey | 


thus be explicitly expressed in terms of known results. We have, on using 
equations (1.6), (1.11), and (2.1), the result 


‘ T : ; aS 
(=) , m1 Pe axe M3). Pils } (5.3) 
Yl w é Ay Py Vo 2X*\ Cn) yo 


hence, on using (1.21), we obtain 
Nu -$,/(Re,,)at{X(1— }.X)}-4 x 


1+ 4(y—1)M3 (ayo (5) 


=U 8116 DAI ag enly-« 


x 


(5.4) 





74 G. POOTS 


In particular, for the case S,, = 1, 


N ee -, pl +dMy—1)M2/a\87-Y297-D as 
Ne fat X(1—4.X)}-4 2? We] (=) . (5.5) 
7=90 


(Re,,) 1+ (y—1)M?2 \a, en 


The ratio Nu), (Re,,) is tabulated as a function of x/x, for a range of Mach 
numbers in Table 6. The quantity 


1+(y 1)? Nu 

1+-4(y—1).Mj5 .(Re,,) 
has negligible variation with Mach number and is given graphically in 
Fig. 4. 


(1V) Re ynolds a nalogy 


teynolds analogy is not applicable to flows with separation. However, 
it is interesting to demonstrate the variation of the inverse of the usual 
Reynolds analogy parameter, i.e. C,Re,,/Nu, given in tabular form in 
Table 6. 


6. Conclusions 


The problem of particular interest is the solution of the compressible 
laminar boundary layer equations for a flat plate in the presence of a 
retarded main stream velocity when the plate is heated. If we assume 
that the Prandtl number is unity and a linear viscosity temperature 
relationship, then using Stewartson’s transformation the resulting equa- 
tions in the incompressible-plane are more amenable to numerical analysis. 
A retarded main stream velocity in the compressible-plane may then be 
obtained on assuming a linear retarded main stream velocity in the in- 
compressible-plane, The essential parameter of the problem is now the 
dimensionless temperature S,, evaluated at the wall. On stipulating a 
particular value for S,. and solving the basic equations (1.8) to (1.10), a 
family of solutions for various free stream Mach numbers M, may be 
derived. The present paper deals with the special case S,, = 1. To examine 
the effect of heat addition to the laminar compressible boundary layer 
when separation occurs it would be desirable to have solutions for further 
values of S,.. This task would entail a great deal of laborious calculation. 
However, the modified Polhausen method developed in section 4 should 
be sufficiently accurate for this purpose. It should be noted that for the 
approximate method the amount of numerical work involved is trivial in 
comparison with that required to obtain the more precise results of sections 
2 and 3 but that its usefulness is restricted to the case o = 1 (when the 
hydrodynamic and thermal boundary layer thicknesses are equal). 

For the boundary condition S,, = 0 equation (1.10) is linear and so can 
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only have as solution S = 0. This will correspond in the compressible- 
plane to a thermally insulated boundary. Equations (1.8), (1.9), and (1.22 
taking S = 0 have been solved by Howarth (4) and Hartree (5). Their 
results will then give data on a family of solutions in the compressible- 
plane for various free stream Mach numbers. Table 5 gives values of z, 
and 7,, (for a range of .M,) corresponding to the case of a thermally insulated 
plate (S, 0) and to the case of a heated plate with heat transfer (S,, = 1). 
It is apparent from this table that heat transfer to the boundary-layer 
moves the point of separation upstream. Furthermore, with regard to 
the family of solutions corresponding to S,, = 1 (see Table 5) it follows 
that an increase in the free stream Mach number leads to an earlier 
separation. Some further facts are revealed from these solutions: 
1+(y—1)Mz2 Nu 
1+ i(y—1) M2 Y(Re,,) 
given graphically as functions of x x, in Fig. 4 are sufficient to describe 
the range of Mach numbers 0 < M, < 6. 

(6) The heat transfer parameter Nu/,(Re,,) for 0 < M, < 6 is virtually 
constant in the region 0 < x < 0-8x, and decreases rapidly in the neigh- 
bourhood of the separation point. 


(a) The dimensionless quantities C,,(Re,,) and 
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APPENDIX 


Relaxation solution of equations (3.5) and (3.6) 


Equations (3.5) and (3.6) subject to the boundary conditions (3.7) constitute a 
two-point boundary value problem. 


Method I for 0 < € < 0-2 


In the usual central difference notation (see Fox (6)), the residual formulae at 
any typical point O, in the range of integration, are 


Ro HU r bhfy) + Fl i— hf) + ( 2 t+-h®gg)'K t (C, thfyC,)% +h*H, (A 1) 
and M, 8,(1 + dhfy) + O5(1 — fy) +(— 2 +h fy)O9 + (Cy + ho Cy )Oo+h* Lo, (A 2) 
where f(n) = + a)O—ady, (A 3) 
dp, 
on a(P- _y), (A 4) 
7 dy ) 


k(n) ~htaS,¥, (A5) 
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Ln) = aS4¥, (A 6) 

H(n) —}aS,', (A7) 

— $89 + 3d%, (A 8) 

AB?+ dhd*, (A9) 

and p is the ‘averaging-operator’ and h is the mesh length. We take as our first 

approximation 0 = S,+S, = 2S, (A 10) 
dd . 1— tf, 

id dn (1 "1 — $4!" 


Point and block relaxation patterns are then constructed for equation (3.6) and 
residuals are calculated at 7 0(0-4)5-2 by means of formulae (A 1), neglecting 
Cw. A second approximation to y% is then obtained using standard relaxation 


and 


(All) 


techniques and the procedure repeated until the first two decimals of x have been 
obtained, Note that at this stage it is sufficient to find ® by using the trapezoidal 
rule. The functions ® and ‘VY are then inserted in equation (3.6) and @ is found by 
relaxation to three decimals. Using this 6 we now calculate new residuals for 
by adding in C'’. The central differences have been obtained near the origin 
using Fox’s procedure (6), and ® calculated accurately. Equation (3.6) is then 
relaxed until we are certain of three decimals in ‘VY. The function @ is then calculated 
accurately using the complete residual formulae (A 2). This repeated relaxation of 
the two equations is carried out until we are sure of three decimals in ¥, ®, and 6. 
The step h = 0-4 is then halved to h — 0-2 and the % and @ equations are relaxed 
so that we are sure of four decimals in the wanted functions. To ensure that the 
fourth decimal is reliable a further run was completed using h 0-1 and working 
to five decimals. The last step is important since we are interested in determiming 
4’ accurately at 7 0. 

The above method was found to be fairly rapid for values of 0-2. However, 
for & 0-2, the initial approximations (A 10) and (A 11) are no longer good, and in 
particular near € 0-6 the coupling between the equations (3.5) and (3.6) is strong. 
Furthermore, the treatment of equation (3.6) for ‘Y becomes difficult since probably 
in this region the non-linear terms are more predominant. It was thus found neces- 
sary to modify the procedure. 

Method II for 0-2 £ Been 

Phe aims of this second method are to weaken the coupling between the equations 
(3.5) and (3.6), and at the same time effectively to linearize equation (3.5). These 
aims may be achieved by using a perturbation approach. We define 

7 .+Y;, (A 12) 
@p ?, Ly, (A 13) 
” 6,44, (A 14) 
1 1 
where ?, |W,dyn and y | Y dy. 
o 0 
The functions ¥;, ®,, and 6; used in the perturbation procedure are defined as 
follows: 


al 
YT, i (ba + Pp’), (A 15) 
dy 


P, ha t db}; ’ (A 16) 


and A, S 4+ SH, (A 17) 
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where dy, 6, and Sy are known functions at the station = €4 and 5), d(¢#)/dn, 
and S\) can be calculated at £ = &, from the series solution given by equations 
(2.16) and (2.17). 

Substituting equations (A 12), (A 13), and (A 14) in (3.5) and (3.6), we arrive at the 
following differential equations in the perturbation functions Y and G: 


ey dY (tb eo (GY aw 
ae | ; } o ah ¢ £ + + - +. — 
dy? | * a1 + a)Oy oa) a - dn * i)} (1 Wan ¥ Lj 


ja¥? 4 a(U%e -U%, e)G4 M,(n) » (Als) 
and 
CG d 
fl | i 
dye t Mi +1 —obadG 


where 


M,( ") 


G iG ' 

Jo, G+4(1+-a)y5——4aGY +N,(n) = 0, (A19) 
Uf] dy 
ay), 
dy? 


dv, (dd, itis er 
ie +a i 44, V+ a(U%e~U%4e)(2+8,), 


(A 20) 


+ {d(1 +-a)D), — ad 4} 


; 120 10; aX. — 
and Ni () - e + {sl + aD, ad 4): , . =" LO, T aS yt Lt 
dy? dyn 2 


dé — > 
HAL bay Ge + ay ¥~ day Y. 
The boundary conditions are 


Y(0) Y(a@) 0, (A 22) 
G(0) (i(x) 0. (A 23) 


Equation (A 18) is still non-linear, but if Y, y, and G are small, the non-linear 
terms will be small and furthermore the coupling between equations (A 18) and 
(A 19) will present little trouble. It was found that the above procedure was ex- 
tremely rapid. In the actual relaxation of (A 18), G was neglected in the first instance 
and Y determined accurately to four decimals using mesh hk = 0-1. The function @ 
was then determined accurately and when inserted into equation (A 18) the new 
Y virtually remained unchanged. Rarely was it found necessary to return to the 
G equation again. It was also found that the amount of numerical work involved in 
solving equation (A 18) could be greatly reduced by the following extrapolation pro- 
cess for Y. By basic iterative procedure we determine at each nodal point four 
successive estimates of Y, ie. ¥(y) for? 1, 2,3, and 4. These are then differenced 
at each value of 7 giving 8,4, Y+)— Y for ¢ 1, 2, and 3. Then we define 
ratios A 5),,/5;-4 for j — 2, 3. Now it was found that A, and A, differed by less 


than } per cent and thus if we define X — 4(A,+A,) then the ‘final’? ¥Y = Y® 8;/(1—A) 
provided A 1. The ‘final’ value of Y extrapolated in this way was found to be 
extremely accurate. Indeed, the fifth decimal was rarely found to be unreliable, 
when verified by a further run, Note that equations (A 18) and (A 19) have been 
arranged so that the relaxation patterns and residual formulae remain unchanged, 
a serious disadvantage of method I. 
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TABLE 4 


Boundary layer characteristics in the incompressible-plane, 
for U, = up(1—4X) and S(X,0) = 1 
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(a) Hartree-Womersley method. 
(b) Approximate solution of the momentum, energy, and thermal-integral equations. 
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TABLE 6 


Boundary layer characteristics for a heated flat plate in the presence of an 
adverse pressure gradient, taking o = land y = 1-4 


For the columns Cy, Rey, Nu/WRey, and C;Rey/Nu the upper figure is that obtained 
using the modified Polhausen method ; the lower figure (in parenthesis) is that obtained 
using the Hartree-Womersley method 
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37. Te = 847 
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SUMMARY 


Solutions of the energy equation of magnetohydrodynamics are obtained for the 
heat-transfer problem corresponding to Hartmann’s velocity profile for forced flow 
between two infinite parallel plates. The semi-infinite plates z +L, x - 0, are 
kept at a constant temperature 7, and the plates z tL, x 0, are kept at a 
different temperature 7, (constant). Solutions are found which are valid for the 
regions « < 0 and «x > O respectively. These are joined smoothly at the plane 
x = 0 by imposing certain continuity conditions. Asymptotic solutions for large 
M 10 are presented, A simplified case valid for large Peclet numbers is worked out 
numerically and the mean mixed temperature and local total Nusselt numbers are 
tabulated and shown graphically. These are compared with the corresponding values 
for the heat-transfer problem in which the magnetie field is absent and the fluid is 
electrically non-conducting. It is found that due to ionic-conductivity the mean 
mixed temperature at any point 1s decreased and consequently the local total Nusselt 
number is inereased., 


1. Introduction 


Tue problem of heat transfer in an incompressible viscous fluid for fully 
developed Couette flow has been studied in a variety of ways by many 
workers, Prins et al. (1), van der Does de Bye and Schenk (2), Schenk and 
Beckers (3), and Dennis and Poots (4). These investigations do not strictly 
apply when the viscous, incompressible, electrically-conducting fluids or 
liquid metals are used as heat-transfer media in the presence of magnetic 
fields. 


In this paper, we consider the simplest problem of this type in which 
an electrically conducting, viscous, incompressible fluid flows between two 
infinite parallel planes under the action of a constant pressure gradient. 
when an external magnetic field of constant strength acts in a direction 
perpendicular to the plates and to the direction of flow. The velocity 
profile for this flow is the fully developed Hartmann’s (5) velocity profile, 
provided that the temperature differences are small enough to permit the 
neglect of convective heat transfer. This is tantamount to the assumption 
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that even with moderate velocities, the buoyancy forces caused by the 
temperature differences are small compared with the inertia and friction 
forces. The infinite parallel plates are kept at a constant temperature 
up to a certain cross-section and after that the temperature is changed to 
a different value. The energy equation of magnetohydrodynamics for this 
problem reduces to a Sturm-Liouville differential equation, Singh (6). 
The mathematical solutions involving two infinite sets of coefficients, valid 
for the two regions, are obtained and the coefficients are determined by 
imposing continuity conditions on the temperature and its first deriva- 
tive at the junction. Asymptotic solutions are given for large values of 
the Hartmann number M (> 10). A simplified solution for large values 
of the Peclet number Pe (> 100) is presented in detail. The first three 
eigenvalues and the associated eigenfunctions are computed for M equal 
to 2, 4, and 10. The mean mixed temperature and the local Nusselt 
number are tabulated and shown graphically. It is found that for a 
constant value of the pressure gradient the local Nusselt number increases 
with M. 


2. Statement of the problem 

Consider a viscous, incompressible, electrically-conducting fluid, moving 
from left to right between two infinite parallel plates under the action of 
a constant pressure gradient in the direction of motion and a constant 
transverse magnetic field H, (Fig. 1). The two infinite parallel plates 
(2 +L, x < 0) are kept at a constant temperature 7, up to a certain 
cross-section « = 0 and afterwards (z = +L,x > 0) their temperature 
is changed to a different value T,. Taking x and z as the longitudinal and 
transverse coordinates with respect to the plates z = +L, the energy 
equation due to Pai (7) simplifies to 


_ oT, eT, &T\ pM ,, Mz sinh M dU) 
We l i = | Fah oe } +- u l z s . = 4 r 
aia? “tamed aut a 7 ts ee? M WV de 


ox" cz" 
(i =], 2), 


sinh M 
M ’ 


with U.= “(cosh M —cosh = A = cosh M — 


4 


where p is the density, c the specific heat, J the mechanical equivalent 
of heat, U,, the mean velocity, k the coefficient of thermal conductivity, 
T, (i = 1,2) the variable fluid temperature in the regions x < 0 and x > 0 
respectively, « the coefficient of viscosity, and M the Hartmann number 
(M = p, LH,(c'p)'), o and p, being the electric conductivity and permea- 


bility respectively. 
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The boundary conditions compatible with the assumption of the steady 
state are 


kel h,T,=0 forx <0; 72. 1,7, —0 forr>Oatz=+L, 
ez oz 

where /, and h, are the thermal transmissivities of the plates. Two 

particular cases which are often met with in practice are (i) h, = hy = 0 

and (ii) h, = h, = 0. The former corresponds to thermally insulated 

plates, the latter to plates maintained at a constant temperature. 


To 











Fic. 1. Flow between two plates. 


The boundary conditions for the problem under consideration are 
T,=T, forx <0,z= +2, (2.3a) 
™%=T, forrz>0,z=+L. (2.3 b) 
The continuity of temperature distribution and its first derivative at the 
junction x = 0, gives 
eT, _ eT 


Cx Ox 


T,=T, and atz=0,-L<z<L. (2.3¢) 


For infinitely large values of x!, the temperature is the particular solu- 
tion of i 


eT, , pM ,, Mz sinhM\  (/dU,\* ad 
kd oa + Ta, | 4, (cosh L yt 2) =% 626 
The solution of (2.1) subject to the conditions (2.3) gives the temperature 


T = T(x,z). The mean mixed temperature and local Nusselt numbers 
are defined by 


L L 
Tig = | Te,2W, de] f Usd, (2.4) 
-—L - ZL 


and Nu = — (=) |m—T,) 
C2 Jem, 


respectively. 
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3. Mathematical solution 
Introducing the dimensionless variables 6, x’, and z’ defined by 
0, = r— T, 
T,—T 
equation (2.1) becomes 
*0, 8 2, ,\)e0; 
in *_Joosh M—cosh u(- 2’— 7 . 


éz2 3n2A\ 


9 \ 
(¢ = 1,2), x = $Pe Lz’, 2= L{-2'—1), (3.1) 


8 


or [sinha (== 1) + A cosh u(=2—1) — A(sinh M) m|, (3.2) 


where 

U Lpc 4 pl 72 
=, . and = a =. 

k = 37Pe @ 4n° A S(T, 


6, and @, are the variable temperatures in the regions x’ < 0 and zx’ > 0 


Pe ——, (3.3, 4, 5) 


T.) 


8 


respectively. The boundary conditions become 
4, 1, for 2’ <0, 2’ = Oand 7, (3.6a) 
6,=0, for 2’ > 0, 2’ = 0and n, (3.6 b) 
fae cH 


6, and ; S. for2z’ =0,0< 2’ <2, (3.6¢) 


C2 Cr 
6, > F(z’, M)+1 as 2’- 0, O : < 7, (3.6d) 


6, > F(z’, M) as x’ > 00, 0: <1, (3.6e) 
where 


= a 
\7i a 


F(z’, M) =| cosh 2a(= = 1] : 84 cosh M{=2 1] 


(cosh 2M +-8A cosh M)+-2'M*(2x—z'){1+-(2A sinh M)/M}}. 
The solution of (3.2) satisfying the conditions (3.6 a, b) is 


° . 
. = v,,(2" sin nz’ - F(z’, ee (3.7) 
qi: . 


3 


n 


2 = ¢ 
*» v,, g(a" )sin nz’ + 2 i.) = pO, Cyc (3.8) 
=i! ° 


3 


n 
Multiplying each side of (3.2) by sin mz’, where m need not be the same 
eigenvalue as n, and integrating with respect to z’ from z’ = 0 to 2’ TT, 


we have 


y » ce. 
‘cosh M—cosh M( 2’—] } ‘sin mz’ dz’ 
| 7 ex 


24*v,,; 

dx’ 

(= 1,2), (3.9) 

with ) x 0 and v,,(0) = 0. (3.10) 


n 
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In order to evaluate the integral on the right-hand side of (3.9), the 
expression {cosh M—cosh M(2z’/7—1)} is expanded in a Fourier cosine 
series in the interval (0,7). The expansion is 


2! 
cosh M—cosh u(= — i) = $b+ 3S 6, cosqz’ (3.11) 
7 a- _ 


om 
=@9 
=, 


: 8M sinh M 
where bg = 2A, Og = — Gant tatty’ 
Therefore 


\eosh M— cosh u(22’ -1)\ 


( Me ‘sin mz’ dz’ = (4+ M sinh ee 
Cr 


M?+-7?n?} dx’ 


S 1672Mmnsinh M dv ing 


rr én), (i =1,2), 
a, [4. M?2+(m— n)*2? || 4.M?2+(m-+-n)?2?| dx’ - mn = ) 


m=1 


(3.12) 


where the summation extends to odd values of m only. Finally equation 
(3.9) becomes 


, dv s 


ni 


Ao 
dx'* ah M?+7 n* 
+ 


x nv 


dr’ : nit 
; 128mnM sinh M dni 
 3A[4M?2+-(m—n)*x* || 4M?2+ (m-- n)*2?| dx’ 


m=1,3 


M sinh alae - 





(ns + 2) = @. 
The solution of (3.13) subject to the boundary condition (3.10) is 


nila’) ni @XPlA, © ~ 
provided that 


— l 28m nAa 


Any+ S . 
Anf(A,m) _ 3A4[4.M2+(m- n)*n? 4ue + (m-—+-n)?x? 


m 


VM sinh M 


mi ' 


(m = n) 0, 


SA 


whuse f,n) = a2? — (4+ M sinh =): - 


327A M?+-7?n? 
The consistency condition for the set of equations (3.15) is 


A(A) aA), = 9, 


2e.2 . 
where —< ot ane . 128yAM M sinh Mo (i 


3A[4.M?+(i—j 2x? if 4 M?*4 -(i+j)?2 2x? | 
and a;, = f(A,t). 


uu 
The determinant A(A) is absolutely convergent and all its diagonal 
elements are quadratics in A; hence A(A) has an infinity of positive and 
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negative roots. The positive roots A,,, (m = 1,3,5,...) are admissible for 
region x’ <0, while the negative roots A,,. (m = 1,3,5,...) are to be 
taken for x’ >0. The two sets of infinite coefficients af (i = 1,2; 
p = 1,3,5,...) corresponding to each A,,, and A,,. are determined in the 
usual way in terms of a‘”"") and a{””) respectively. Now we write 
2S Q ; ye — 
0, > exp(A,,. ¢ UZ ny (2’ )+ roe (z ’, M)+ l (x : 0), (3.17) 
m=1 
25 ¢ : 
0, = - EXP(Am22")Zna(2’) + +5 F(z’, M) (x’ >0), (3.18) 
7 


m=1 
where Zinil2’ = 5 aj sin pe’ (6 = 1,2). (3.19) 
1 


The two sets of coefficients ates and a\””) are determined by the conditions 
(3.6c) at 2’ = 0, namely, 


9 >! ae 
=: = S ay") sin pz’ i . aA?’ ,M)+1 
hase A 


= 
7 


9 all a 
_ b 3 > al"®) sin pz’ + ~, F(z’, M), (3.20) 
1 


m=1 p 
x x 
and Y YA ay"” sin pz’ :. s Ana tp sin pe’. (3.21) 
m=1 p=1 m=1p 


Multiplying both sides of (3.20) and (3.21) by sin mz’ and integrating with 
respect to z’ from 0 to 7, we obtain 


> a (p = 1,3, 5....), (3.2% 
1 


3 


Pp m 
* (m1) (m2) 
> Ami 4) -> Jn ; 


The two equations (3.22) and (3.23) dalieaine av”) and a’), Because of 
the slow convergence of the left-hand side of (3.22) we make use of a 
variational principle. The eigenvalues and the associated eigenfunctions 


can also be calculated for other forms of the boundary conditions (3.6 a) 
and (3.6b). 


4. Asymptotic solution for large values of MV 
For large values of M (M > 10), both cosh M and sinh M tend to $e”, 

and hence 
8 | M sinh M) 8 { M? 


+ Pa... _\ 
372 Al mn?+ M2 ~~ 5a (M—1)(n?n?+- M®)| (4.1) 
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and the expression under the summation sign of (3.13) becomes 





eh M*mn ni ees (sn) 
(M—1)[4.M?+-(m—n)*x?|[4M?2+(m+n)*x®}] | M 
and may be neglected. This corresponds to the case of a magnetic field 
ef high intensity. The equation (3.13) then reduces to 


od*v,; 8 f M* \dv,, .¢ , 
2— ni a a — 1,2), (4.2 
‘ “de? Sai) (M —1)(7?n?+-M?*)} dz’ went > 
The solution of (4.2) subject to the conditions (3.6) is 
v, (2) = a, ,exp(A,;2’) (t= 1,2 (4.3) 


where 7? BAnif 4 = ial —n? — 0. (4.4) 


"i 3a2 || + i 1yntn® M)| 
Hence the complete solutions for the two regions are 


9 = 


ad hy . *\ 2 ~ Q Vf! i P 
= Pi a,,, exp(A,,, 2’ )sin nz’ + rE 1(2,.M)+1 (mn 
n=1 


qo = 


6, = ss > a,9eXP(A,,9 2’ )sin ne’ +@ Re, M) (x’ > 0). 


n=1 
The continuity conditions give 


2 
anit ” = Ano, 


Ani @n1 = Ana @ns- (4.8) 

Hence msneeee ee. 
ence a.3 = n(l Ace)’ and a,,= wh(l—AnAa) 

Table | gives the first three eigenvalues and the coefficients in the case 
M — 10, corresponding to Peclet numbers 1, 10, 100, 1,000, and o. In 
‘Table 2 are given the eigenvalues corresponding to Pe = 0 when M = 15, 
20, 50, 100, and o. The local Nusselt numbers for these values of M at 
x’ x are given below. 





yi : 20 5° 100 @ 


| Nu]x ‘275 ° 2°358 2°420 2°443 2°467 


TABLE 1A 
Eigenvalues and coefficients for x’ < 0, M = 10 





I 10 100 1,000 D 





3°3256 168°46 16,520 16 x 10° 
79072 185°73 15,914 16x 10° 
12°5800 218°36 15,568 15 x 108 
06688 070384 00004 
02961 0°08 39 O°0013 
0° 1869 0°09g02 0°0022 
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TABLE 1B 
Eigenvalues and coefficients for x’ > 0, M = 10 





I | fe) 100 1,000 D 





1°6709 3°2954 3°3588 3°3612 
6°3189 "733 31°400 31°460 
11'032 “550 87°410 89°640 
I°3212 *Qg615 1°9996 2°0 

0° 3706 “5828 0°6654 0°6666 
O°2131 »* 3098 0°3978 0°4000 




















TABLE 2 
Pe L 





is 50 100 





34041 ¢ 3°6270 3°6641 
7 


31°69 32°6 32°59 


59°47 go"d4 g2°§18 














5. Solution for large values of Peclet number 
For large Peclet numbers (Pe > 100), a? in equation (3.13) is of the 
order 10-° and can be neglected everywhere except in the immediate 
neighbourhood of the points « = 0, z = +L, where the plate tempera- 
tures change abruptly, giving rise to infinite temperature gradients. 
Near these points it will not be permissible to neglect the derivatives 
2*v,,,/Cx"* in (3.13). The analysis of this section breaks down near x = 0, 
+L. With this proviso, the diagonal elements of the determinant 
(3.16), instead of being quadratic, are linear factors in A and hence the 
infinite determinant has only an infinitude of negative roots valid for the 
region x’ > 0. This excludes the possibility of heat diffusing into the 
region x’ < 0. The temperature distribution in the region 2’ < 0 is inde- 
pendent of the x coordinate, This idealization makes the solution valid, 
to a good approximation, only for high Prandtl number fluids and when 
the velocity U), is sufficiently large. The expression for the temperature 
distribution in the region x’ < 0 is 


6, = SA’, M)+1. 


The boundary conditions for the region 2’ > 0 become 
6, 0, for 2’ > 0, 2’ = 0 and 7, 


Q 


6, gehil? M) ‘1, fora’ =090 


i , ' 
0, —> “AM asz’>oo (0 <2’ <7). 


~ 
= 
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The solution of (3.2) satisfying the conditions (5.2 a, c) is 


@ 


5) ee 
6, = = > exp(A,,, 2’) pa a,” sin pz’ +7 a 5 F(z’, M), 


m=1 


where A,, are the roots of the determinant 


where 
128ijAM sinh M 


“ 3A[4M?+-(i—j)?x*|[4M?-+ (i+ j)?x?]’ (@ #4), 4 = f(A,t). (5.5) 





The coefficients a!” can all be expressed in terms of aj” and, using the 
condition (5.2 b), we obtain 


(m) . s 5.6 
> am = = (5.6) 


m- 


These determine a” but because of the slow convergence of the right- 


m 


hand side of (5.6), we make use of the variational principle and get 


D 


met 3 x ) 
in) ¥ . (m) gin) re 
> — > Te ie (5.7) 
1 


m=1 p- 
yp 
The mean mixed temperature and the local Nusselt numbers are evaluated 
from 
16M?cosh M < as” 


7A Z omG..# ) > P(r’ pt +4 M?)’ 


m 


and Nu = — z2 exp(A,, x’) 2 pay”. 


In order to calculate the eigenvalues approximately we study the 
determinant formed by retaining the first four columns and the first four 
rows of (5.4) and evaluate the first three roots A,, (m = 1,3,5) for M = 2, 
4, 10. The corresponding set of coefficients a)” (p = 1, 1,3, 5) are calculated 
by the use of the variational principle (5.7). Tables 3, 4, and 5 give the 
first three eigenvalues and the set of constants, the mean mixed tem- 
perature and the local Nusselt numbers. 
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TABLE 3 


Eigenvalues and coefficients for M = 2,« = 0 





m I 3 





Am 2°9105 








0022 
0°072 
0° 489 











Nusselt number for M = 2,« = 0 





Nu 





4°69 
“86 
‘70 
18 
03 


-~_— =— — = BH NH Ww 











TABLE 4 


Eigenvalues and coefficients for M = 








90°96 





(5) 
ay 





0056 o'o2!I 
0579 o'1og 
o°0s4 0450 














Nusselt number for M = 4,« = 0 





6 Nu 


m 





~O4 4°75 
792 
2°79 
“29 
14 
09 
06 


wv NN NKR 
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TABLE 5 
Eigenvalues and coefficients for M = 10, « = 0 





I 3 5 





— 89°44 


(3) (5) 
a ay 








ool4 0°008 
0638 0033 
o'o1g 0421 














Mean mixed temperature and Nusselt number for M = 10, « = 0 





x’ 6 Nu 





4°96 
4°20 
309 
2°54 
2°37 
2°29 
2°254 
2251 
2°251 
2°251 











6. Comparison and discussion of the results 

The eigenvalues and the coefficients of the cup-mixing mean temperature 
found for various values of M are given in Table 4. The mean mixed 
temperature and the local Nusselt numbers are plotted against x’ for these 
values of M in Figs. 2 and 3. It can be seen that the local Nusselt num- 
bers go on increasing with M for the same values of x’. Fig. 4 is the plot 
of [Nu],, against M. For M tending to infinity the Nusselt number tends 
to 2-467. Hence it can be concluded that the heat transfer coefficients of 
the liquid metals are increased because of the presence of the magnetic 
field and the thermal entry lengths are considerably decreased. 


TABLE 6 
Comparison of eigenvalues and mean temperature coefficients for 
various values of M 





Prins Present authors 


M=o M =2 M=4 M = 10 








1°885 1°9403 2°0438 2°2407 
21°43 21°6133 21°0 20°938 
62°31 57°54 60°64 59°627 
O'ol4 0909 0893 0° 869 
C053 0069 0066 0078 
Orols oo16 0026 o’o3I 
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T - = rv 
. . °-8 
0-4 o6 L, 
2. Mean mixed temperature against x’. 


























‘1G. 3. Nusselt number against 2’. 
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20 3or 40 
M— 


Fic. 4. |Nuj, against M. 
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DISPLACEMENT FUNCTIONS AND LINEAR 
TRANSFORMS APPLIED TO DIFFUSION 
THROUGH POROUS ELASTIC MEDIA 
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SUMMARY 


It is shown that the stresses and the pore pressure in porous elastic media through 
which a liquid is diffusing can be expressed in terms of two displacement functions. 
These functions are particularly useful in problems relating to a semi-infinite body 
or infinite layer when the stresses or displacements are prescribed on the surface. 
Problems of plane or axially symmetric strain are closely related when expressed in 
terms of these functions and the introduction of suitable repeated transforms 
(Fourier or Hankel followed by a Laplace transform) permits a parallel development 
of solutions for these two types of strain. 

The process is time-dependent since it involves diffusion. After an infinite time 
the pore pressure reduces to zero and the governing equations for the stresses and 
displacements in the medium reduce to the ordinary equations of infinitesimal 
elasticity. The limiting forms of the functions used in the paper can then be related 
to known stress functions. 

The Fourier and Hankel transforms can be regarded as specializations of the 
double Fourier transform, and it is shown in conclusion how the plane strain analysis 
ean be generalized into a three-dimensional treatment. 


1. Introduction 

THE typical problem with which this paper is concerned can be stated 
thus: a liquid (usually water) diffuses through the pores of a clay medium 
under the action of stresses applied at the surface, and we seek to determine 


the displacements, the stresses, and the pore water pressure in the medium. 
Problems of this type are embraced in a theory which workers in the field 
of soil mechanics term ‘the theory of consolidation’. We assume at the 


outset that the medium is elastic but the stresses and the excess pore 
pressure (i.e. the excess of the pore-water pressure over the hydrostatic 
pressure at any point—in the text, we often drop the qualification excess) 
are time-dependent since diffusion is a time-dependent process. A peculiar 
feature of these problems is that we cannot isolate the evaluation of the 
pore pressure from the evaluation of the elastic deformation: the two must 
be derived together. It will be seen below that many of our equations are 
formally identical with the equations for the determination of stress in 
an elastic solid due to a non-uniform temperature distribution and the 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960] 
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displacement functions which we shall presently introduce can be used in 
the solution of thermal stress problems. The analogy is not however com- 
plete, and, in general, thermal stress problems are much simpler since the 
temperature distribution is either given or determinable without reference 
to the stresses. 

We shall limit ourselves here to the explanation of methods which 
facilitate the task of solution (solutions of some outstanding problems in 
soil mechanics are given elsewhere). 

The equations which govern the diffusion of pore water in clay were 
given earlier by Biot (1) and the present work is based on Biot’s analysis 
which we shall quote with a minimum of comment. An adequate explana- 
tion of the physical content of the mathematical analysis will be found in 
most textbooks on soil mechanics, e.g. Terzaghi (2). 

Problems in soil mechanics frequently relate to the semi-infinite body 
or the infinite layer, and in succeeding sections of this paper we shall 
examine separately plane strain and radially symmetric strain. It will be 


shown at a later stage that the two analyses can be united. An observation 
of minor interest in the present context is that the time-dependent equa- 


tions of our analysis reduce to the classical equations of infinitesimal 
elasticity for t-+ «. It may therefore be anticipated that the displace- 
ment functions used in the paper are related to known stress functions of 
the theory of elasticity and this point is investigated briefly. 

Biot (3, 4) has recently recast his earlier work (cited above) in a more 
general form. Most of these generalizations envisage a more complex 
physical situation than we shall be concerned with; but one of his conclusions 
may be quoted here since it is pertinent to our work: he remarks that 
‘from a mathematical standpoint what distinguishes the consolidation 
problem from an elasticity problem is the addition of the function ¢ 
satisfying a heat conduction type equation’ (4). We have reached a 
similar conclusion by a different route. 


2. The governing equations 

2.1. Notation 

In the main, we have followed standard notation but we have introduced 
a few changes which conform with present usage in mathematics and soil 
mechanics. An unusual feature is that soil mechanics workers prefer to 
regard compressive stresses as positive and we have followed them here. 
The principal symbols follow. 

t Time coordinate. 
Z, 2, p; U, W, U, Coordinates and displacements in the horizontal, 

vertical, and radial directions. 
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Volumetric compression per unit volume. 

Total compressive stresses. 

Excess pore-water pressure. 

Effective stresses acting on the clay skeleton. 
Shearing stresses. 

Shear modulus and Poisson’s ratio for the medium. 
An auxiliary elastic constant. 


Coefficient of permeability in Darcy’s law. 
The coefficient of consolidation. 
L< e? 


ép* ; p Cp oz? 
2.2. Plane strain in a fully saturated medium 
Biot'’s equations (1.3) (1) written in the above notation are 


> 1 C 
V2u—(2n—1)— ——— = 0, 
Cx G Cx 
Vw (27 - = 0, 


ce 


et’ 


cV%e 
take compressive strains as positive and replace 1/(1—2v) by 


From these we can show by differentiating the first two of equations 


(1) that V2(e-4 2(' ne) = 0, (2) 


and, if S is any harmonie function of x and z, we can write 


Fal 
— ne ° 


+-Uy+-Ug, 


If we now write 


uw Ww, ¢ UW's + Ws, 


the differential equations (1) are satisfied if the following equations are 
satisfied: 
, OF : l OW 
V7u,+ ’ Vw, + >= —(¢ +4 (4a) 


Cx ox C2 


+ To avoid a possible misunderstanding, we remark that equations (1) to (4) and (7) on 
pp. 290, 291, of (2) do not appear to be consonant with Biot’s analysis: we make no use of 


them in the present paper. 
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CU, , CW, 
—*+—#=0; (4b) 
Cx Cz 
OUg , OW, 
and V2w, —— + ——* as 0. (4¢) 
a 2 
It is readily verified that a solution of equations (4) is 
¢Y oX 
V2X, u= ——, wy 
Cx 
os . 
WwW, = z——S, 
oz 
oY ey 
a... v2 
Oz 
where X is an arbitrary differentiable function. If we now write 
E=X+Y, 


then . co 
CE es 


ce od -» 


Y 
Y 
—S, 


at 
2a(- — iN*8), 
C2 


w 


and the governing equations for F and S are (cf. the third of equations 
(1) and the definition of S) 
. ,CE , o 
cV4E = V?—, VS = 0. (7) 
ot 
We can now write down expressions for the stresses in terms of EF and S 
but we defer this for the moment. 


2.3. Axially symmetric strain in a fully saturated medium 
The displacement equations for cylindrical strain are easily derived 
from Biot’s analysis: 
1 ce 1 @oa 
V2 ——Ju,—(2n—1)— —- — = 0, 
p* ép G Op 


ée 1 éa 
V2w—(2n—1)——-— — = 0 
C-- Gh 


cV2e = - 
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102 


where 


and 


we can show, as in section 2.2, that 
V2(o-+ 2(rne) 0, (9) 


The functions £ and S can be introduced by a procedure similar to that 


used in section 2.2, and we find that 


ck os 
( p ; cp 
. OF 
cViE = V?—., 
ct 
Ex} ressions for the stresses in terms of Eand 8S 

The strain of the medium is related to the effective stresses by the 

ordinary equations of infinitesimal elasticity (the pore pressure causes no 


strain), e.g 
cu 


2G] ——+(n iy, 


For comparison, we write down in parallel columns the expressions for the 
displacements and total stresses when the strain is plane or axially 
symmetric. 


Plane strain Axially symmetric strain 


rap Dy rar! Ck es 


u=: — z P = - Ba» 
cr cr cp cp 
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‘\p--S+S 


dz? ° Gz 


Pe i ie 


: E —-—+—, 
pop 


zéS . as 
pcp cz 


Oe Al “27 a2Q 
o. é ie _ 8h aS 


20 Oxéz 2G pez pez 


It can be seen that the shear modulus G appears in these expressions only 
as a multiplier and the elastic ratio » appears only in the expression for o. 
This is convenient in applications since the elastic constants are introduced 
at a fairly late stage of the analysis. It can also be seen by inspection 
that—with the exception of ogg—the expressions for cylindrical strain are 
immediately obtained from the corresponding expressions for plane strain 
if x is replaced by p everywhere and the appropriate form of V? is used. 


2.5. The use of non-dimensional quantities 

It is convenient to anticipate here the application we wish to make of 
the formulae of the preceding section. The typical problem may be stated 
thus: a uniform pressure f is applied to a strip of width 26, or to a circular 
area whose radius is 6, on the boundary of the medium. The complete 
specification of the problem will then require the intervention of four 
constants (, f, b, and ¢ (the coefficient of consolidation), together with the 
dimensionless ratio ». The five constants can be reduced to two if we 
choose 6 as unit of length, 6?/c as unit of time, and f as unit of stress. 
The components of displacements and stress are then replaced by dimen- 
sionless quantities, but the expressions tabulated in subsection 2.4 are 
unchanged in form, save that the coefficient ¢ must be replaced by unity. 
To demonstrate this point completely would require a complete change 
of notation and we shall indicate only the first few steps. We write 


x be’, p Os, s=be’, (= tH, 

using a prime for the new non-dimensional coordinates. Similarly we 
write E = bE’, S = bS’. 

The differential equations for FE and S are transformed to 

4 , 20K’ Qo 
Vik’ = V?—_, V2S’ = 0. 
ot 

The displacements are replaced by u’ = u/b, ete., and the stresses by 
o,,/f, ete., though it is often convenient to retain the non-dimensional 
forms o,,/2G. We can retain the existing notation and drop the primes 
provided it is understood that non-dimensional quantities are used. We 


shall presume this understanding in the following sections. 
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3. Fourier and Hankel transforms of FE and S 


The results of the preceding sections enable us to treat the theory of 
plane and radially symmetric strain in parallel. We remark that the 
introduction of the functions F and S reduces the number of dependent 
variables from three (e.g. u, w, a) to two, but the order of the governing 
equations has been increased. The reduction may prove to be an illusory 
advantage unless we can at the same time eliminate one or more of the 
independent variables by means of suitable transforms. The method is 
familiar and may be stated briefly. We assume that EF and S can be 
represented by definite integrals and that these integrals exist and can 
be differentiated the required number of times; in specific problems it is 
usually possible to verify existence and differentiability a posteriori if a 
solution can be found. The differential equations governing E and S are 
transformed into simpler equations governing two new functions, say F, 
and S., the boundary values of E, and S, being determined by an inversion, 
The more common practice of defining transformed functions E, and S, 
and later transforming back to E and S can be justified a priori but is 
less convenient here since it requires the transformation of the displace- 
ments and stresses—or at least of those displacements and stresses which 
appear in the boundary conditions. A minor but not unimportant point 
is that the introduction of new symbols for the transformed stresses is 
an additional burden on an already overburdened notation. 

The Fourier cosine transform is appropriate to problems of plane strain 
in which the loading is symmetric about the z-axis. We represent the 
stress functions FE and S by the definite integrals 

x 


B(z,2,t) = — | cos(2€)E(E, z, t) dé, 


. 


2 
. 


| cos(x€)S.(E, 2, t) dé, 
0 
and similar integrals for the displacements and stresses can be obtained 
by differentiating under the integral sign. 
The zero order Hankel transform enables us to eliminate the coordinate 
p in problems of cylindrical strain. We define E and S by 


E(p,z,t) = [ EJy(p&)E,(E,2,t) 4€, 
0 


x 


S(p,2,t) = [ EIy(pe)Sy(E, 2, t) aE. 


0 
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If, for example, we form the derivative @S/ép and assume that we may 

differentiate under the integral sign, we find that 
es a : 
a ES, (pE) S, dé. 
cp - 
i) 

If we differentiate again and make use of the differential equation for Jp, 
we find 


V?S(p,2z,t) = | é(08)| —# |S dg 


= 0. 


This equation is true for 0 < p < and, if we multiply by pJ)(fp) dp and 
integrate from 0 to «, we have by Hankel’s theorem (see, e.g., (5), 
Theorem 135) as 

(5 : #)S\620 =0 (0<é<~). (13) 
The expressions for the stresses and displacements may be tabulated very 
compactly if we use the notation 


K(%,4) = =cos(ee) or Ed(pé), 
7 


p 


K’(f,¢) = —Zesin(xg) or —E2Ny(06) 
p 7 

and accept the convention that the first kernel on the right of these 
equations is to be used in plane strain and the second in axially symmetric 
problems. We use a single symbol u to signify the horizontal or radial 
displacement. 


= | (7. §)[—Beat28.n] a6 
0 





@ 


o 


-(x ,\aS 
0 


In problems where surface tractions are specified, expressions for the 
components of total stress are required and these are derivable from the 
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above formulae for the displacements and the pore water pressure. For 
example _ 
Cun 2 7 : Ow Ad es Y 

et k( 16) | eh + — 2628] dé, 
A i p 
0 


CZ 


i 2) 
— = _(x ,\[eE eS 
z,%e — [ K'(*,2)| 29° ag. 
2G 2G J p Cz Cz 
0 
The differential equations for E,, and S,,, are 


o =) e(e.2,0) _ 0, 


(14) 
€) 81620 = 0. 


In these equations, and in the table above, the quantities E and S should 
bear the appropriate suffix c, or h, as we have indicated in the expression 
for wu. 


4. Repeated transforms of E and S 

The partial differential equations (14) can be transformed into ordinary 
equations by means of a Laplace transform. Remembering that the initial 
state of the medium is determined by the consideration that the pore water 
diffuses with a finite velocity, we know that the volumetric strain is zero 
initially. Since e— VE. 


we must have : 0 fort = 0; 


or, equivalently, ——€7)E., = for t = 0. (15) 


We now drop the suffix c, h except when necessary to avoid ambiguity. 
If we write ' ytio 
_ | EE,2,p)dp (y > 0) (16) 


youn 


E(é,z,t) 


PA a) 


then, using (15), 


/ 


sh cee 
oe I co eas 
se ses oe FAL 
a 2 > 

ct\ cz* 


and the first of (14) becomes 


(= ¢2 
ae ¢ 
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This equation is true for all positive values of t and we shall assume 
E(é,z,t) =0 (t <0). 
If we now multiply the equation by e~?dt and integrate from —« to , 


we have by an obvious modification of Fourier’s double integral theorem 
[see (6) equation 7.1 (2)] 


2 2 = 
(i e eer 0, (17) 
dz? 12 


the left-hand side being analytic in the half-plane re(p’) > y. Similarly, 
the second of (14) becomes 
(va —-#)5 - 0, (18) 
dz* 
y+ia 
| ertS(é,z, p) dp. (19) 
yrix 
It is unnecessary to rewrite the table of section 3 in terms of double 


where S(€,z, t) ; 
aT 


integrals since we can in most problems perform the integration with 
respect to p before using these expressions. 


5. Bilateral displacement functions for infinite layers 

The functions F and S of section 2 are convenient for the solution of 
problems relating to the semi-infinite body z > 0 when the loading is 
applied at z — 0. They may also be used in problems of the infinite layer 
0 <z<—h, but it is sometimes preferable to modify these functions so 
that the boundary conditions at the upper and lower surfaces may be 
simply incorporated. Since three conditions must be satisfied at each 
surface, the boundary conditions form a sixth-order matrix and it is 
expedient to simplify the algebra as much as possible at the outset of the 
work (a boundary condition on the pore pressure o must be satisfied in 
addition to the usual two boundary conditions of elastic theory). 

We shall illustrate the use of bilateral displacement functions by con- 
sidering briefly their application in plane strain problems. Retaining the 
notation of section 2.2, we write 

u=u'+uw’, 
E = E'+4 
It is easy to verify that, if z, 
, ok’ . os’ 


‘= — 


‘ = ‘ , > 
ox ex 
C E »i C S - 


Ox 1 ox 
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satisfy equations (1) if o is given by (3) and if FE and S are governed by 
(7). 

These functions may then be expressed by a repeated transform and 
the solution can in general be written as a linear function of the terms 


exp(—z&), exp[—2(€*+-p)'], exp(—z,&), exp[—z,(€*+-p)*] 


of which the first pair reduce to unity for z = 0 and the second pair when 
z h. 


6. The limiting forms of F and S (t-- ~) 

The steady state solution of the diffusion equations which we have been 
examining is identical with the solution of the corresponding elastic 
equations. The excess pore pressure reduces to zero in the steady state; 
equations (1) and (8) then become the familiar elastic equations for plane 
and radially symmetric strain. We may then expect the limiting forms 
of E and S to be closely related to known stress functions of the theory 
of elasticity. We shall show how F is related to (i) Love's biharmonic 
stress function y which is applicable in problems of symmetrical strain in 
a solid of revolution (7), (ii) the harmonic stress function of Boussinesq (8) 
for the solution of problems relating to the semi-infinite body (9). Love 
(10) notes that equivalent formulae were given by Hertz (11). 

Denoting the limiting forms of F and S by E,, and S,, and writing 


(V4L = 0), 


el 


C2 


we have, from (3) and (6), S, AVAL 


if we ignore linear terms. The displacements u, and w become 


VL, 


Cx 
,eL 


Cen 


and v 
If we express ¢ in terms of Love's stress function y, we have 


?, ? 
; 1—2 2©X_ 
2G Cz 


Since ¢ is harmonic, we may write 
I ov Ci, 
4: [-" @, 
aq * 


where w is an harmonic function which can be determined by comparing 
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the expressions for the displacements; we ignore here an arbitrary bi- 
harmonic function of e alone since the function may be absorbed in y. 
We find finally that 


where 


and the expression in square brackets is harmonic. 

A comparison with Boussinesq’s stress function ¢ (we follow Green and 
Zerna’s presentation here) is most easily obtained from an examination 
of the plane strain case. The displacement is given by 

eL 


CXC 


u 


D eae 
a on 12 . V?L, 
Cr 


and the expression for w is unchanged in form. If we write 


ks 
» ~ 1 2 _ 

L = 2L,-z (V2Z, = 0) 
it is easily verified that 
€ 4 
29 » 1 
7 (27 1)z Dad 
CXLCZ 
and the shearing stress is zero on the plane z - 0. For comparison we 
write side by side our expressions for u and w in terms of L, and 


Roussinesq’s expressions in terms of ¢: 


anc. we may write 


7. Three-dimensional displacement functions in consolidation 

problems 

The comparison with the stress functions of Boussinesq, Hertz, and 
Love, suggests that we may be able to carry the preceding analysis a stage 
further by formulating techniques for the treatment of more general types 
of surface loading on a semi-infinite body z > 0. We shall outline briefly 
a quite general method of attack by making use of a special form of 
Fourier’s double integral. The use of this method in elasticity problems 
is well-known (12) and it will appear that the solution of a consolidation 





110 J. McNAMEE AND R. E. GIBSON 


problem is in general feasible if the corresponding elastic problem can be 
solved. 

It is easily verified that the plane strain equations for EF and S and the 
equations which express the stresses and displacements in terms of these 
functions can be generalized into three dimensions without change of 
form, the plane operator V? being replaced by 


. - 
cy” 


For definiteness we shall consider a problem which frequently arises in 
applications of the theory; an area D of the plane z = 0 is subjected to 
a normal stress f(x,y), the shearing stress being zero all over this plane. 
In addition, a boundary value of the pore pressure o, or its derivative 
6a éz, or a linear combination of both will be prescribed; it is sufficiently 
illustrative here to take this condition to be a = 0 on z = 0. We adopt 
as possible solutions for E and S 


"+4 Aye ay +PM)eilar+Bv+rl dadBdp, (20) 


| | | Byetiarst-ve-et dadBdp, (21) 
v . *. . 
- 2 » Br 


where y x24 B2 


77 


and Br denotes the Bromwich—-Wagner contour. The conditions of zero 
shear and zero pore pressure on z = 0 yield two relations between the 
functions A,, A,, B,; then from the formulae of section 2 the boundary 
value of the normal stress o,, can be expressed by 
| | AA, etor+Pu+m dadBdp (22) 
“os Br 
with A 1+-ns—(l+s)? (8 


Inverting this relation we may write 


“ 
/ 


TS feed ical +B dtdé. 
“ 4n*p J J 2G 


D 


For example, if f is constant and D is a rectangle of sides 2A, 2, the origin 
being at the centre of the rectangle, we find 


yAA, ae 
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The functions F and S in this example are found to be defined by 


24 p L - f sinaAsin Bu, 
j = 


fe-=71+#)t (1 +4 g)he-*7 Jeux +Bur+7t dad Bd p, 
pAy'*aB 


—o-o Br 


@ 


of [ [ sin aA Sin Be oars 8y)-2+0" dadBdp. 
f 27rt : Ay*aB 


The solution which we have just sketched is purely operative and we have 
not considered at any stage the validity of our operations. It is clear, 
however, that the existence and differentiability of the Fourier integrals 
in (20) and (21) depend on the behaviour of the Bromwich contour integrals 
qua functions of a and 8. Equation (22) shows that the only singularities 
in the p-plane are poles and branch-points and it is usually possible to 
effect the integration with respect to p. 

The detailed exploration of these techniques is considered in a following 
paper. 


. 
—-o-—-o Br 


8. Conclusion 


The solution of specific problems of the theory we have been considering 
will lead in general to algebraic work of considerable complexity and very 
few solutions are available. The use of displacement functions followed 
by repeated transforms permits an easy and natural incorporation of the 
boundary and initial conditions. The success of the method depends 
mainly on the possibility of evaluating certain integrals in the complex 
plane, the integrands of which possess, in general, a finite number of poles 
and branch-points. 
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SUMMARY 
This paper gives a method for the calculation of the influence functions arising 
in Multhopp’s subsonic lifting surface theory. It is eminently suitable for pro- 
gramming on automatic computers but double length working would be necessary 
if large values of the arguments and high order functions were required. 


1. Introduction 
CERTAIN aerodynamic theories initiated by Multhopp (1, 2) involve several 
functions of two variables whose numerical values are required in practical 
applications. A few of these functions have been tabulated by Curtis (3), 
to supplement charts given in (1), Figs. (1)-(6). Desk machines have now 
been replaced by automatic computers in the rest of the calculations 
involved in the theory and a suitable method for forming each value of 
the functions directly in the machine is needed. The purpose of this paper 
is to meet this requirement for a reasonably large number of Multhopp’s 
influence functions, and for a wide range of the independent variables. 

The functions considered depend essentially on the quantities 
* [eos p9+-cos( p+ 1)6][ X — $(1—cos 4)| dé 

{| X —4(1—cos @)|?+ Y?}3 


, l 
S(p) | 


T(p) | | cos pO +cos( p+ 1)0}{[ X —4(1—cos 6) |? +2} dé 


0 





and are given by 
(X,Y) = 14+ S(0), j(X,¥) = 4S(1) 
k(X, Y) = 8(2), (X,Y) = S(3), m( X,Y) 
X 
ii( X,Y) | (Xo, ¥) dX, = X—}+7(0) 
X 
(X.Y) = | j(X,, ¥) dX, = 14+47(1) 


x 





[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960) 
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x 
kk(X,¥) = | k(X,,¥)dX,_ = 7(2) 


x 
W(X,¥)= | WX, ¥)dX, = 7(3) 

2. The method of computation 
The functions (2) are not linearly independent; each can be expressed 
as a linear combination of any three others. If i, j, ii are known, suitable 
expressions for deriving successively the remaining six functions are as 

follows: 

12k = 8i—3j-+32ii—8.M(4i-+-j —8ii)— 128 i+ 128Y2, (3) 
= 8i—j—4k, (4) 
321 = —24i+3j—4k + 96ii—48.M(2i-+ 2k—jj)—96Lj, (5) 
96kk = 3j—4k—8l, (6) 
15m = —12i+4k—I-+ 12jj—12M(j+ 41—8kk)—192Lk, (7) 
48I1 = 4k—I—3m, (8) 
where M = 2X—-1, L = X?4+Y?—X-+3. (9) 


The proof of these formulae is given in Appendix A. 

Each of the influence functions can be expressed as a linear combination 
of three complete elliptic integrals. If the following three elliptic integrals 
are chosen, relatively simple expressions result for i, j, ii. 


Let 
in 
2 [ ai [1—(P—Q)? Jcos*a d ae 
7} {4PQcos*a+[(P+Q)*— 1 ]sin?a}! 
dn 
“ dx 


= | RP Coo [POF Tp 


$77 


[ =. PP Q?) sin? dx 

r [4PQ cos*a+ (P+ Q)*sin?a]{4PQ cos*a+ [(P + Q)?— 1 Jsin®a}t | 

(10) 

where P = [(X—1)?+Y?]#, @ = [X?+ Y?}!; then, as is proved in Appen- 
dix B: 

i= 14+A+4(2X—1)B+C, (12) 

j = 8(1—X)A—16Y2C, (13) 

ii = (X—})+- (4X4 })A4[2¥?+ (2X —1)(X—})]B+(¥24+X—})C. (14) 

Bartky (4) gives a method of calculating complete elliptic integrals 


5092 .40 I 
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which is eminently suitable for an automatic computer. This iterative 


method is an extension of Landen’s classical transformation. Bartky’s 
result is that if 
in 
Ili) F (a,m,cos*6+-b,r,n,;sin?6) dé 
1 —$ —_______ 


(m, cos?6+-r, n, sin?6)(m? cos*6-+-n? sin?6)* 





0 
= 4(m,;+n,), Nis, = (m,n,)* 


1.+7.0. 
h(a; } b,), ” cain octets 
TM 
nN. 
b+ (r, 1 
4m; 4, 


then I (i) I(i+1) 





If we now put g? = r;m;n, and assume r;, m,, n; are all positive, we 


l m,n; 
Vi+1 aC aie 
: [™s nj " 


To adapt this result for the present purpose of computing A, B, C, we 
write 


can deduce that 


(16) 


Further, 9;., = m;., if q; = m,. 


my = 2[ PQ}, ny = [(P+Q)*—1}* 
Ag 1—(P—Q), by a= @, Co = 0 
dg = P—Q, % = P+Q 
and let bes 4(m,-4 


\ (17) 


= (n,m,)! 


4(n;a;-+-b;) 


3(Me, td. 
z qi . ; 


4i 





The iteration proceeds until m, 


q, to the accuracy of the calcula- 
tion. Then 


B= 


1 C me Orth. 


m, m,, 


(19) 


for a,,, d,,9 and ¢,., c 


pag When m, = n, = q,- 


3. Comments on the method 


Some of the coefficients in the relations used increase as |X| and Y 
increase resulting in loss of precision by cancellation. Higher-order func- 
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tions (i.e. those with larger p) are affected more than lower-order ones. 
The cancellation can be overcome by using multiple-length arithmetic 
but a large range of applications is adequately covered by single-length 
working. When |X) and ¥ are as large as 3, for example, only three deci- 
mals are lost in forming II. 

The convergence of the iterative process (18) is rapid and three or four 
cycles of calculation are usually sufficient to give six decimal places 


correctly. 

Errors will be introduced when P or Q are small, that is, very close to 
Y =0,X =0orY =0, X = 1. For then m, is small. Theoretically the 
method also breaks down when Y = 0 and X lies between 0 and 1, for 
then n, = 0 and so the final m, = 0. It has been found, however, that 
correct values may be obtained if the rounding error in the formation of 
square roots is such that the calculated root is always strictly positive, 
for then n; # 0 and the process converges. To ensure this and also to 
avoid trouble with the scaling of numbers the computation of the iterative 
process is best arranged as follows. The multiplication routine should 
produce a double-length product; the division routine should work with 
a double-length dividend and a single-length divisor to give a single-length 
quotient; the square root routine should produce a single-length root from 
a double-length argument. Then if the precision of quantities b; and d; 
is double-length and of the remaining quantities single-length, no further 
thought need be given to the size of the numbers occurring in the iterative 
process. Moreover, it will be found that difficulties do not arise when Y 
is small. For instance with Y = 0 and X = 1-000004, we obtain i with an 
error of only 4x 10-°. 
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APPENDIX A 
THE LINEAR DEPENDENCE RELATIONS 


Equations (4), (6), (8) may be proved by integration by parts as follows. Putting 
R(@) = X—4(1—cos@), W(0) = (R*+ Y?*)t, (20) 
we have e . e 
cosnOW(0) d0 sin nO 1 [ sinn@sin OR(B) dé 


rol | 
; O)) +3 | n WO) 
0 0 





(21) 


n 
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where the first term on the right is zero if n is a positive integer. Further we have 
identically 

lsinn@sin@  Lsin(n+1)@sin@  cos(n—1)6—cos(n+1)0 | cosn8—cos(n+ 2)6 
2 n ‘s n+l 4n 4(n-+1) 


1)6@+cosn@ cosnO+cos(n+1)@ cos(n+ 1)0+cos(n+ 2)6 








cos(n 


4n : 4n(n +1) 4(n-+-1) 





The sum of the two equations derived from (21) with n land n = 2 now gives 
equation (4). Similarly we deduce equation (6) from n 2 and n 3, and equation 


(8) from n 3 and n 4. 

To prove the remaining formula we make use of the identity 
R(O) 
W(6)” 


Substituting from (20) for R(@) on the left-hand side and for [W(@)}? on the right- 


cos mA R(b)W(8) cos mO{ W(8) |? 


hand side we deduce 
4M cosmé W(8)— }[cos(m + 1)8+ cos(am— 1)6]W(8) 
RC RO 
Leosmé r a { }M[cos(m 1)8+ cos(m 1)6} Ti . 
R(8) 
20j——. 
lire) 


The integral with respect to 6 from 0 to 7 of the sum of the two equations derived 


99 
22 


A; [cos(m + 2)8 + cos(m 


from (22) with m O and m 1 is 
4M(ii— X + }) > }[ii-—- X+}34+}(jj—1)) L(i— 1) +}3M(i—1+ }j)+ &(4I+k). 
Substituting in this equation for jj from (4) we obtain (3). Similarly (5) may be 
obtained from m l and m 2 by substituting for kk from (6), and (7) from 
mi 2 and m 3 by substituting for Il from (8). 


APPENDIX B 
THE TRANSFORMATION TO ELLIPTIC INTEGRALS 
We use a transformation similar to that given by Curtis (5) whose equation (3) 
is equivalent to our (29) below. His method does not appear to be suitable for 
deriving all the equations we need. Our transformation is 
(P—Q)—(P+Q)cosa 
(P+Q)—(P Q)cos a’ 
where a ranges from 7 to 0 as @ ranges from 0 to 7. This relation gives 
{PQ sin ada : 
Y)—(P— Q)cos «]*’ 
4PQ sin? 
[(P~Q)—(P— Q)cosa]?’ 
PQ{4PQ costa + [((P + Q)— 1]sin®a} 
{(P+ Q)—(P— Q)cos «]* 
To simplify our subsequent equations we put 
d P-—-Q, f= P+Q \ 
Gla) 4PQcos*a+[(P+ Q)?— 1)sin®a | 


H (x) G(a)+sin®’a f?—d* cos*a 


cos@ (23) 


sin 6 dé 


sin?@ 





(we)}2 





MULTHOPP’S INFLUENCE FUNCTIONS 
Equations (24), (25), (26) yield immediately 
do _ _ 2W(6) 
or. ~  [G(a)]# 
‘7 
| A rf 
2a} W(8) “ted [G(a ii | [G(a)}* 


and thence 


Equation (28) can also be used to prove that 


0 T 
] ( co0 a0 a “d—feosa 2dxa 4 [ 9tfoses 2 dg 
27) Wie) f—deosa (Ga)? 22) f+ dcosp (GB) 


0 0 
where B a7—a. Taking the mean of the last two expressions we derive 


1 [cos dé l | 2fd(1—cos*a) da 
2a | ~Wi(8) 22 | f?—d*® costa [G(a)}4 
0 


A final application of equation (28) gives 


7 0 wT 
lL { sin?6 dé | | 8PQ sin*a dx I F sPQsin?B dB 


a} W(8) a | (f—deosa)*® [G(a)}A f +dcosB)? (G(B)|! 


Taking the mean of the last two expressions, we derive 


7 T 

l | sin*9 dé l [ SPQ sinta (f* +d" costa) f : x ; 

a) W(@) 7. (f?—d? cos*a)? [G(x)}? 
0 in 


Now 


Aes COs asin a} 
dx H(x) 


1 
[ G(x) 8H (a) }? 


{2 cos*a sin’ad*G(a) + (2 costa — 1 )d*G(a)H (a) — 


cos?a sin?’ad?(1—d?)H(a)}. (34) 
The expression in the curly brackets on the right of (34) is a polynomial in cos*a 
of the third degree which can be expressed in terms of the independent quantities 
[ H(x)]* cos*a, [H(a)]*, H(a)sin*a, sin’a( f? +d? cosa) in the form 
[1 — (1 —d?*)cos*a){ H(a)]* — sin?a f2d?H (a) — ( f?—d*)sin*a( f? +d? cos*a). 


Using this expression, integrating (34) with respect to a from 0 to 47 and substituting 
from equations (29), (31), (33) we obtain 


7 
1f R(6) dé fas 
- cos@ Wid) A. (35) 
0 
We are now in a position to prove the relations connecting i, j, ii with A, B, C, 
namely (12), (13), (14). Equations (29), (31), (35) give (12). Next, we have 
; [sin 0W(4)] 
1@° 


1 
: wa! }M + (X*2+ Y?—X)cos 8 + } Mcos*@ + 3 cos*6] 
Integrating this result with respect to 6 from 0 to 7 we obtain an expression for 
* cos*? ; , ‘ . : 
we which, where used in conjunction with (31), (35) gives (13) and also with 
{ 


(29), (31), (35) gives (14). 
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SUMMARY 

An approximate method for obtaining a solution of the simplest form of the 
transport equation is described, The approximation is of the ‘diffusion’ type, but 
it incorporates fitting certain parameters by least squares, and so, it is hoped, gives 
better solutions than the straightforward approximation of this type. 

A numerical solution is exhibited for a simple case, and this shows features 
characteristic of ‘diffusion’ solutions, but gives a good qualitative and quantitative 
solution. 

The method can be extended to obtain an approximate solution of the transport 
equation in other coordinate systems. 


1. Introduction 
Tue simplest case of the transport equation, the one we will designate 
‘one-dimensional’, is 


1 
potinjeS a, Pal) | ile Po) dy (1.1) 
C2 r=0 4 
Here 7 is the intensity of radiation at a point in a semi-infinite isotropic 
distribution of scattering particles in a given direction. The point is 
defined by the single space coordinate z (a pure number measuring the 
depth below the surface in multiples of the mean free path of the particle 
distribution), and the direction in which the intensity is measured by @, 
the angle between the given direction and the normal to the surface z = 0. 
As usual, » = cos 6. 

The optical properties of the scattering material are given by the ‘albedo’ 
for scattering w (< 1), which measures the fraction of radiation scattered, 
and the scattering function 


f(O) = s , P.(cos®), (1.2) 
r=0 


which measures the density of radiation scattered in a direction at an 
angle © from incident unit intensity of radiation. 

+ Present address: Mathematics Department, Sir John Cass Coilege, University of 
London. 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 1, 1960] 
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The usual boundary conditions are that the intensity be specified over 
the forward hemisphere at z = 0, and over the backward hemisphere at 
some other point z = Zp. 

Using the Wiener-Hopf technique (1), an analytical solution is possible 
when the latter boundary condition is given at infinity. In general, how- 
ever, the complexity of the analytic solutions is so great that it is 
necessary to find numerical solutions also, To obtain any reasonable 
degree of accuracy in solution, the time and space required is such that 
it is only really feasible using a modern high-speed computer. 

For this reason there has always been considerable interest in approxi- 
mate methods. The most important and successful is that of Chandra- 
sekhar (see 2), which replaces the integral on the right-hand side by 
a quadrature formula. Solving a system of linear differential equations, 
and using the properties of the Legendre function, a solution accurate to 
any desired order is obtained. As the corresponding treatment for the 
equation when expressed in spherical polar coordinates is not so con- 
venient, a method was sought that could be applied to more general forms 
of the transport equation, 

Another method involves solving what is, in fact, a different physical 
problem by the use of the diffusion equation. This method gives solutions 
which are least accurate near the boundaries—usually the places of 
greatest interest. 

The method proposed here is based largely on the latter method, but 
should in general give a better answer. 


2. The integral equations 


We shall consider, for simplicity, the case when only the first two terms 


appear in the summation on the right-hand side. Bearing in mind the 
normalization condition implicit in the definition of the scattering function 
this then becomes 
: ; 
jen | i(z, po’) dy’ H Bu | i(z, mw’ dy’ > (2.1) 
1 1 


where £ is a positive constant, This represents an ovoid scattering pat 


tern, that is, one scattering more radiation forward, and less backward, 
than in the isotropic case, 


We shall take as boundary conditions: 
a(O, 2) 


and i(Z», pe) = O 
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Because these conditions are only specified on half boundaries it is 


convenient to extend the notation. We put 
0(z, ) I(z, p) if 
J (z, pt) if 


and write further 
1 


p(z) | {1 (z, pe’) 4-J(z, p’)} dp’ 


| {I(z,y')—J (z,p’) iu’ dy’. 


0 


and q(z) 


Then equation (1.1) becomes 


O<z< 2 and O<u: 


el 
we +1 = dol plz) + Bya()), 
and Ps -« hen| — p(z) + Bugq(z)]. 


ez 
The corresponding boundary conditions are 
1(0, ht) 1 
J(Z,) = 0 for all p. 
2.2) and (2.3) as if the right-hand sides of the equations were 


Integrating (2.2 
known functions, we obtain integral equations for /(z, 4) and J(z, 4), viz. 


and 


Hz, p) = 8h | eben) + Bugle) dt (2.6) 


0 
Zo 


and J (z, p) = | ec OE p(t) — Byq(t)} dt. 
“HY 


2 


1 
| p! np -a\p dys E,(a), 


t) 


whence, by the definition of p(z) and q(z), we find that 


Now 


Zo 
plz) = E,(z)+40 | [ p(thE,{ t—z\} 4 Bsgn(z—t)q(thk,f t 
6 


z\}| dt, 


Ey(z) + 4en [ [sgn(z—t)p(t)E,f t-—2|} + Bg(t)Egf\t—2 }] dt, 


0 


q{z) 
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a pair of simultaneous integral equations for p(z) and q(z). Ifw = 1, then, 
from (2.7) 


Zo 


co-=-0046 p(t) —Pug(t)} dt 


| e™-2)4{ p(2g—U) —Bug(zo—u)} du, 


and therefore 


I (z, 1) +-J (zg—z2, p) = e-F/#+ | ef 4[ { p(t) +-p(zo—t)} + 
) 
+Bu{a(t)—q(zo—t)}] dt. 
It is easy to see that this integral equation is satisfied identically by 
T(z, )+-J(z)—z, uw) = 1. (2.10) 


3. Differential relations 


If we integrate (2.2) and (2.3) between » = 0 and p = | and subtract, 
we obtain 


(3.1) 


and if we multiply (2.2) and (2.3) by », integrate and add, we obtain 


= 


dz 


dr ( po 


\ == 0, (3.2) 


where by analogy we have put 


r(z) = | Ue, H')+J(z,p’)jm? dp’. 
0 
We may obviously continue this process to obtain differential equations 
for the higher moment functions. There will, however, always be one less 
equation than there are moment functions. 


4. The approximate solution 

We shall now assume that : 
r(z) == Ap(z), (4.1) 
and substitute this into (3.2). There is no obvious immediate physical 
interpretation of this assumption, though it implies a diffusion-like solu- 
tion. The justification lies in its utility. 

Eliminating q(z) between (3.1), (3.2), and (4.1), we find 

d*p 
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where y = 1—}8u. If thenmw + 1, 
p(z) = Ae*+ Be-™, (4.3) 


and q(z) = Ax + Be az__ Aer}, (4.4) 
¥ 


where Avt—y(l—aw) = 0. (4.5) 

Now whatever values are given to A, B, and a, if (4.3) and (4.4) are 
substituted into (2.6) and (2.7) we have values for I(z, 4) and J(z, 1) which 
satisfy the boundary conditions. 

Thus we will choose A, B, and « to give a good solution of the pair of 
integral equations (2.8) and (2.9) in the sense of Gauss. 

Let f(z) and g(z) be the functions obtained by substituting (4.3) and 
(4.4) into the right-hand sides of (2.8) and (2.9). Now 

Ze Zo 
| sgn(z—t)q(t)E,{|z—t|} dt = A sgn(z—t) aP py z—t\} dt 
J YJ dt 
0 0 

&o 
| — | pOVEy(it—2)} t+ 2ple)—pled)Eleo—2) —P()E2)], 

"L 
after an rie by parts. Hence 


of =A, 


pe) (0) Bale) +55 = Pleo) Haleo—2) —= ple) + 


2 
+ jo( ote | jm (t)E,{\z —t |} dt, (4.6) 
and similarly 


q(z) ! - Te no) E;(2) 3M plea) —2) + a1 


w 
2(1—a) 


Ze 
; r atl bal 
< [q(%) Ee(zo—2)-4 alO)Ex(2)—2402)] + 55 (1+ 2) [ g(t)E,{\2—t\} dt. 
(l—) vt 
0 


(4.7) 


Zz 
The integrals r,(:,2) = [ eWE(t) dt (see 2 and 3) 
0 


can be determined, and so 
Zo 


eM“E{\t—z|} dt = * tanh Ages 4 1 (1 -+a)2}— Ey{(1 —a)(z9—2)}] + 
a a 


~25) —~ Ble). (4.8) 
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Thus if (4.3) and (4.4) are substituted into (4.6) and (4.7) the term in 
p(z) which for large z is of order e® is 


aa _ BA tanh tx BA) oa: me , Ba -ar)| 
Ac| (1 -} Aw} 1 t 


x y) ye 


/ 


te 


Similarly the term in q(z) of order e™ is 


Ao, _ BU _ tanh hax) 


e xz 

x | ; x 
Thus if « is chosen to be a root of the equation 

Bil tanh~! 
on ee ( wy x) _ 
| 2 | x | 

the terms in p(z) and f(z) which are dominant for large z will be the same, 
as will also those in g(z) and g(z). Since the left-hand side of the above 
equation is an even function of a, the coefficients of e~* will then also 


0, (4.9) 


agree. It can soon be verified that, for small values of 1—za, there is a 
single pair of symmetrically placed real zeros in the strip —1 < rea < 1. 
The appropriate value of A can be determined from equation (4.5). The 
roots of (4.9) for some values of w and 8 are tabulated below. 


TABLE 
Root & of the equation 


_ Bd w)\{,__tanh~ta) 
ae A x | 





0°93750 0°90625 





259555 0°405259 "490590 
"291471 0°407425 "493074 
"293045 0° 409549 “495540 
0°294017 o'411663 "498005 
2960177 07413705 *SOO451 
0°297729 0°41 5557 
0°209274 0°417938 
> 300810 0° 420009 





* 302338 0°422070 











°84375 031250 78125 





oOO1s14g oO O6507 4 712986 
*O2Z10O14 0°671624 "715968 
YOH2 3505 o'674559 Y7ISO3S 

0°677479 *721887 

0°680385 *724524 


0°68 3277 


i) 27746 


o' 6861 $4 30054 


7 
7 

o'HSgo1d "733547 
7 


640076 o691 868 36426 
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The coefficients A and B are now determined from the condition that 
the mean square deviation of (z) from p(z), and of ¢(z) from q(z) shall be 
a minimum. 


Let H(A, B) = [ [{ p(z) — p(z)}? +{¢(z) —q(z)}*| dz (4.10) 
0 


where p(z) and q(z) are given by (4.3) and (4.4). Then 


CH _ 9, oH 9 
cA cB 


are a pair of simultaneous equations to determine A and B. 

When these values of «a, A, and B are substituted into (4.3) and (4.4) 
we shall have an optimum solution of the integral equations (2.8) and 
(2.9). If these functions in their turn are substituted into (2.6) and (2.7) 
we have finally 


- 4) — etl — tel 4 _(p Pm) Bf, BU) |) 
I(z,p) ! j r —() : 1) 4 es —(14 - '4)) | + 


J 4 ( POR lens B (: , BU- ™ u)e] (4.11) 
1+ ap x l—ap a 
and 


“ A f, B(l—w) \,., B {, Bili—m)| i 
J (z, hor] - 1+ plex c | Bon eet 1 
Gm) =% r= : x a T 1+ opl x i° 


A | B {, Bil—a)) | - 
; | oto nt oe pase ele-20H, (4,12 
™ l ye x ? whee I 14 xu i‘ € ( ) 


If the region extends to infinity (zg = 0), the appropriate solution 
would require that A = 0, since the intensity is bounded. 


5. The case of no absorption 
When w = 1, equations (4.2), (3.2), and (4.1) give 
p(z) = Az+B 


and q(z) = .% A. 


Y 
Then, if we define p(z) and g(z) as before, 


p(z) = Az+ B+ E,(z)+4A 1 +9 s(2)— Bale —2)} 29 A) 7 


—4B{E,(z)+ E,(zy—z)}, (5.3) 
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qe) = —y4{+¥) 4 
:. ¥ 
+E,(z)+4$A { l +" \E Qe) + E,(zo—2)} +2 a) se 
/ P 
- } BLE; (z)— E; (2 —2)}. (5.4) 
The terms in j(z) and p(z) which do not decay exponentially agree 
whatever the value of A. The constant in ¢(z) is the same as that in q(z) if 
A=}. (5.5) 
(It can easily be verified from (4.5) and (4.9) that A 4 asw—1. This 
is in fact the well-known Milne-Eddington approximation.) 
As before, A and B are determined so as to minimize 


H(A, B) = { [fp z)}* + {q(z)—q(z)}*] dz. 


Now £,(z)-+E,(z)—z) and £,(z)+ E4(z,—z) are even functions of z— 4zp, 
and E,(z)—£;(z,)—z) is an odd function of z—4z 9, and so 

Zo 

| {E,(z)+ E,(z9>—z) }{ E,(z) — £3 (zp—z)} dz 

0 


Zo 


] {E;(z)— E3(z—z)}{ E,(z) + E,(zy—z)} dz 


0 


and | {E3(z)+ E3(z)} dz = | { 


0 


Hence the equation 


reduces to }(Az,)+B 


which is equivalent to equation (2.10). 


6. Numerical solution 

A numerical solution, using the above approximation, has been com- 
puted for the case when w = 1, 8 = 0-1, and z) = 2-0, as a solution of 
the complete equation with these parametric values already existed but 
is unpublished. It had been computed by solving the integral equations 
(2.6) and (2.7) by a Picard iteration method. 

The comparison of the two solutions is shown in Fig. 1; the approxi- 
mate solution is drawn as a full line, and the solution of the full equations 
as a dotted line. 

It will be seen that the greatest discrepancy is always near to each 
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Absorbing wall-—— 


Approximate solution 


Numerical solution 








die 4 4 4 
o4 08 12 i) 2-0 


> 
~ 





Fic. 1. Solution of the one-dimensional transport equation, with linear 
scattering, when there is no absorption, and 8 = 0-1. 
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boundary, a fault which, as pointed out above, is characteristic of ‘dif- 


fusion’ approximations. The approximations get progressively worse as 
the direction in which the intensity is measured moves from the forward 
direction towards the side. Even at its worst, however, the approximation 
gives both a good qualitative and quantitative description of the solution. 

The author wishes to thank Miss K. M. Stocks of the Armament Re- 
search and Development Establishment, for her help with the computa- 
tions. 
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